{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cox Processes (w/ Pyro/GPyTorch Low-Level Interface)\n",
    "\n",
    "## Introduction\n",
    "\n",
    "This is a more in-depth example that shows off the power of the Pyro/GPyTorch low-level interface.\n",
    "\n",
    "A Cox process is an inhomogeneous Poisson process, where the arrival rate of events (the intensity function) varies with time. We will use a GP to model this latent intensity function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import torch\n",
    "import gpytorch\n",
    "import pyro\n",
    "import tqdm\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create sample training set\n",
    "\n",
    "Here we're going to construct an artificial training set that follows a non-homogeneous intensity function.\n",
    "\n",
    "#### Intensity function\n",
    "\n",
    "We're going to assume points arrive following a sinusoidal itensity function.\n",
    "This will give us points in time where there are lots of arrivals, and points in time when there are few."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "intensity_max = 50\n",
    "true_intensity_function = lambda times: torch.cos(times * 2 * math.pi).add(1).mul(intensity_max / 2.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Determine how many arrivals there are"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of sampled arrivals: 48\n"
     ]
    }
   ],
   "source": [
    "max_time = 2\n",
    "\n",
    "times = torch.linspace(0, max_time, 128)\n",
    "num_samples = int(pyro.distributions.Poisson(true_intensity_function(times).mean() * max_time).sample().item())\n",
    "print(f\"Number of sampled arrivals: {num_samples}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Determine when the arrivals occur\n",
    "\n",
    "We'll sample the arrivals using rejection sampling. See https://en.wikipedia.org/wiki/Poisson_point_process#Inhomogeneous_case for more details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def log_prob_accept(val):\n",
    "    intensities = true_intensity_function(val)\n",
    "    res = torch.log(intensities / (true_intensity_function(times).mean() * max_time))\n",
    "    return res\n",
    "\n",
    "arrival_times = pyro.distributions.Rejector(\n",
    "    propose=pyro.distributions.Uniform(times.min(), times.max()),\n",
    "    log_prob_accept=log_prob_accept,\n",
    "    log_scale=0.\n",
    ")(torch.Size([num_samples]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The result\n",
    "\n",
    "Here's a plot of the intensity function and the sampled arrivals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEFCAYAAADkP4z+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO2dd3wc1bXHv7NFvXfbsq3ibuOCx70bTDE4dAhJnEAgJCEkz5DQXkgIPJI4JA8ImOQl1NBrwBiIK9jG3ePeZUm2ZMuqVm+rLfP+2F0hC5VdaXdny/1+Pv54Nbsz97d3Zs7eOffccyRVVREIBAJB6KDTWoBAIBAIfIsw/AKBQBBiCMMvEAgEIYYw/AKBQBBiCMMvEAgEIYZBawGu8NBDD4nQI4FAIOgDy5cvlzpvCwjDD/DYY4/1ed+KigrS0tI8qMYzCF3uIXS5h9DlHsGo69FHH+1yu3D1CAQCQYghDL9AIBCEGMLwCwQCQYghDL9AIBCEGMLwCwQCQYghDL9AIBCEGD4L55RleRbwIaAC84GrgQogXlGUFb7SIRAIBKGOL0f884EBiqIMAFKBZEVRXgcSZVme5q1GVVXl9xvO8M+vTrP7dA0iDbWgv1isNr44UcmvVx7l/g8P8z/rivn7plOU17dqLU0QBFQ0mFh1sMyrbfhkxC/LchpwLXCnLMt3AQuAY463jwKLgZ09HaOioqJPbZfVt7H6eA2rj9cAcPGgaJbNGUR2ckSfjudJamtrtZbQJUJX9/z7UBVv7Kmkssl8wfa1ebWsP1rK328YppGyb+IP/dUVQlfXtFltvLu/iteUclotKsNirSRHG72iyyeGX1GUCmCKLMtjsbt7NgM1jrdbgYzejtHXlWsRsWYeXthIcSN8dricvSVN3P7eSZYtzOWuOVl9OqYn8ceVgiB0dUeNuZbKJjNZyVFcP3EA6fER1NTWsb+8jQUjUjTX1xl/0+NE6LqQovPN/Pjd/ZyqagZgdm4SMfGJpCVFeUWXT1M2KIpyRJbll7GP+KMcm2OB895qMy7SyOLRSaSlpfGLBbn89YsC3lbO8r/r82k1W/nFwlxvNS0IQu5fNIzZw5KZnZuEJNlToFRU6Ll93oU35p/WnCQnJYqbJg/SQqYggDhV1cT3X91LRYOJ7JQofrt4JDNzk73apk98/LIsd0wS1AY8AYx3/D0GWO0LHQlRRh69ehR/vn4cep1EXETApCoSaESr2cr/fHac6qY2AAx6HXOGJbcb/a44XtbAy9uK+O2qY2zKq/KVVEEAUljVxNJX91DRYGLK0AQ+vGuq140++G5y90ZZlrfJsvxLYJOiKFuBVlmWbwdqFUXZ7CMdACwZn8GnP5vObTOH+rJZQYChqir//fFR3th1lmXvH3J5v1EZsdw9LxubCvd+cIgT5Y1eVCkIZOpbLOgliWnZifzze5OIDvfNYNRXPv73gfc7bXvCF213R05KdPvr8noTZquNzMRIDRUJ/I0XthTx2eFyosP1/PrKkW7t+4sFORSdb+azw+X85M39fHL3dGLFE6agExMHx/PxT6cRbtATFab3Wbshv4Dr4Nk6rvn7Dpa9f4g2i01rOQI/oaCyiWe/LADgqRsvYmR6jFv7S5LEH68dw7iBcZyra+XPa096Q6YgQHG6DgESo8J8avRBGH6GJkcRadRzqKSef245rbUcgR9gs6k88slRzFaVmy4eyPwRKX06TrhRz/LrxmDUS7y7p4Rdp2t630kQ9FQ1mrj82W38ac1JrDZt1hWFvOGPjzSy/LoxALyw5TSldWIRTqjz9u6z7C2uIzU2jAcuG96vYw1Pi+HHc7KYMjSBtNhwDykUBDJPbyigvtVCfmUjuu5jBLxKyBt+gGnZSVwxNo1Ws42/rBOP5KFOcU0Lep3Eo1eNIi7S2O/j/XRuNq/dNpms5KjePywIao6cq+fDfecw6CT++4oRPUaHeRNh+B3cv2g44QYdnx4qZ0+xf64sFPiGh68YwRfLZrFotGcWzRj0OnQdhnZaPd4LtEVVVX7/nzxUFZZOG0x2hwATXyMMv4PMxEjumGUP7/zj6jyR0yfEyYj3fEqP/IpG7nx9H/+7Pt/jxxb4P+uOVbKnuJakaCN3z8vWVIsw/B340ewsFo1OZZlYzRuS/GPzKdYfr/Daj36rxcZX+ed5Y+cZyutNXmlD4J+oqsrzmwoBuGdejkdciP1BGP4ORIXpWfHtCczuZWWmIPgoOt/MX78s5OfvHORsTYtX2hg3MI7Lx6Rhstj4++ZTXmlD4J80mawMTYoiPS6cGy8eqLUcYfh7wmIVcf2hwktbi7DaVK6dOIDBSd6bhP2vhblIEnywt4SqRjHqDxViIgw8e8t4Pr9nBuFG38bsd4Uw/F1QWtfKsvcO8bN3DmotReADqpva+PhAKQB3zsryalu5qdEsHJmK2ary9u6zXm1L4H/E+CglQ28Iw98F4QYdG/Mq2ZhXxbHSBq3lCLzMO8pZTBYb84Ynk5vq/UiL22cMAeCt3WdpNVu93p5AW55en8/2wmq/ChgRhr8LkqLDuPFiezrdN3ed0ViNwJu0WWy8tcs+8r7dR0n75KEJjB0QS6PJyqGSep+0KdCGwqom/u+r0/zkrf00mvznR14Y/m74ztRMAFYdKqO+xdzLpwWBymeHy6hsbGNkegzTsxN90qYkSfzh2jFsvHc2U7J806ZAG95xuPOuvijDr5L0CcPfDTkp0czMSaLVbOOj/aVayxF4iUtGpvLgZcP5+YIcn0ZyjcqIJTkmzGftCXxPS5u13XZ8Z0qmxmouRBj+HnCO+t/afdav/HMCzxEXaeSHs4Z6bJWuu1isNo6cE+6eYOSzw2XUt1qYkBnH2IFxWsu5AGH4e2DBiBTS48Iprm7meJkophFsaP1j3tJm5dK/buWWF3dT09zW+w6CgEFV1fa5I38b7YMw/D1i0Ot48vqxrF82i9EDYrWWI/AgrWYrS/62g6fX52u2XiMyTM+w1BjMVpVPDpRpokHgHY6UNnCktIGESCNXjk3XWs43EIa/F6ZnJzEoQVTmCjbWHK3gZEUTWwurMei1uw1ummxfxfnB3hLNn0AEniMtNpyfL8jhjllD/WLBVmf8Z5rZz1FVldoWM4lRYkIuGPhg7zkAbtJ4+fyCEakkRRvJq2jiUEk94zPjNdUj8AxpseHcMz9HaxndIkb8LpBX3sgVz23np28d0FqKwAMUVzez63QNkUYdV43L0FRLmEHHNeMHAF//GAkE3kYYfhfITIykosHEvjN1FJ1v1lqOoJ98esjuT180Oo0YP4itvn6S/alj9dFyUfc5CHhy7Un++dVpapv9d/2PMPwuEBWm5/Ix9nC/lQdETH8go6oqnx4qB+Dq8dqO9p2MSI9hRHoMMeEGztZ6JzOowDfUNpt5bUcxT2/Ix+THP+LC8LvINRPsj+MfHyjFJiooBSwnyhspqGwiIcrIzJwkreW08/LSSWxYNoscDasyCfrPf46UY7aqzMxJIj3Of2ssa/+cGyBMy0okIy6cktpWDpbUM3GwmIQLRIanxfDqDy6mosGEUcNons6kikLsQcFnDjfiEsdA0V/xnyvfz9HppHZ3z+oj5RqrEfQVvU5iRk5S+xOcv1HRYKKgsklrGYI+UNFgQimuJcyg49KRqVrL6RFh+N3gCsdCjB2nqjVWIugL/h4nv+ZoOXP/9yv+vO6k1lIEfWDd0QpUFWbnJvlF0EBP+Lc6P2NiZjwvf38SU0VGxYDkT2tOcup8M79YkON3uVMAJg9JQAK25J+nvsWseV1WgXusPmr3BFw5zv9W6nZGjPjdQKeTmJWb7Fe+YYFr2Gwqnx8pZ2NeFVY/HfmnxIQzJSsRs1Vl48kqreUI3OR70wazeFw6C0f4t5sHhOHvM40mi4juCSAOnaunvN5ERlw4F/nhaN/JZY4soeuOVmisROAul49J5+mbLvJ7Nw8Iw98nfvfpcWY8uZmDonpSwLDWYUgXjU7zad59d7l0lH20uDn/PC1t/lOxSRBcCMPfB4x6iTaLjbXHxKgsEFBVlXWOc3XZGG3y7rtKRnwE4wfF0Wq2sbXgvNZyBC5wvrGNhz8+wsa8wHHP+dTwy7I8SpblzxyvfynL8lJZlu/xpQZP4ByVfXGiUmMlAlc4WdFEUXULiVFGJg9J0FpOrywanYZBJ3FKpAcJCL7Mq+Tf+0oDqj63zwy/LMvhwGVAtCzLs4FkRVFeBxJlWZ7mKx2eYPKQBOIjDZyqaqawSsRc+zvOJ7NLRqWi1/mvm8fJLfIgtj0wlx/NztJaisAFvjhhH+kv9PPY/Y74chbiduBF4HpgMXDMsf2o4++dPe1cUdF3t0ptbW2f9+2OaYNjWJtXyyfKKb5zcd/cB97Q5QmCTdf8IWHoZg9kVHpkv66j7vBWf1U09G//YDuP3qYvulrNNrbk2w3/+GQpYK4vnxh+WZYvBb5SFKVZlmWAFKDG8XYr0Gu2rLS0/vlm+7t/Z66aqLI2r5adJS0su6Lvx/a0Lk8RTLrS0mCcl1Oje6O/rDaVsvrWfhUCCqbz6Avc1bXheCUmi8q4gXGMzRnkJVWe7y9fjfh/BKQ7jP5EYB6w3vFeLBBws1izhyVj1EvsP1NHTXObKNAi8Chna1q48Z+7iArTs2HZLL+ORAplnPN8l4xK0ViJe/jE8CuKcovztSzLG4FfA1cC7wFjgNW+0OFJYsINPH3TRYzKiBVG3495/LPjRBr1/GDGENICKBHawPgIdJJESW0rJyuaGJEeo7UkQSesNpUvHf79S0b55xNMd2gSzqkoylagVZbl24FaRVE2a6GjvywancbgRFGP119pbrPy/t5zvLStiACY070AnU5i3ohkgIAKEwwlTBYbN08exLzhyYxIC6x02j5fYqYoynzH/0/4um1voqqqeBz3M7YXVtNmsTEhM46UmMAZ7TtZMCKVf+8r5cu8Su6ak6W1HEEnosL0LLskV2sZfUIs4Oonr+0o5qoV29lWKDJ2+hvOkfL8EYHlf3UyMzepfR6puqlNazmCIEIY/n5S1dhGfmUTm8TjuF+hqipf5tkn3hYEQNKsrogJNzA1KxGbCl/lB1z8Q1BT2WDilW1FAbuORxj+fjJ3uH00uemkuDH9iaOlDVQ2tJEeF86ojMCdGHU+rew6XdPLJwW+ZNPJKpavOcmTawKzdoL/p5HzcyZmxhEXYeD0+WaKq5sZkhSltSQB9hsTYP7wlICee7lqXAaTBicwdkCs1lIEHdjsGOg5B36Bhhjx9xODXsesXHv0xWYx6vcbpmYlceuUTK4YG1hhdp1JjgnjokFx6AItLCmIMVu/TqA3d3iyxmr6hjD8HsAZdrdJFM/wG+ShCfzu6lHMzA3MG7MrRP0H/2DfmToaTVZyU6PJDNBwbmH4PcCcYXbjsvNUjcihLvA4p6qa+M5LCj/41x6tpQigPZBjXoCO9kH4+D1CSkw4912Sy8iM2IDI/hjsvL+nBINeYuHIVOKDoG5takw4+8/WAdDQaiE2ACo8BTPOCKtA9e+DGPF7jB/PzWb+iBTCDKJLtURVVZ7fVMhDHx3lbE2L1nI8QkyEgYmZ8VhtKjtPifUiWmK22shMjCQ1Jiwgajt0hxg6CIKKgsomSutMJEeHMTojeCJhZuUmsae4li0F1Vw6OrAnrAMZo17H326dgM2mBvSEuxieepCVB0q59/1DVDWatJYSsjgfw2cPSwroG7Mzsx3zSKIco38Q6NeWMPwe5JODZXx+uJzthWKxjVZscRj+OcMC1//aFeMG2teLFFe3UFwtSjJqgaqqbC+sptUc+AEcwvB7kFk5SQBsE6MyTWiz2FCK7dWKZjrORbCg10nMcHynrQXCz68FhVXN3PavvVz53HZUNbBDa4WP34PMzHXcmIXVIlunBuw7U0ur2caI9BiSY4KvRsJ3pmaycFQqc4JobUIg4XSzXTwkPuDvbWH4PciItBiSo8MorzdRWNVMbmpg5egOdFRgytAEJg6O11qKV5ieHVxPMYGGMwNvMCwKFIbfg+h0EjNyEvn0UDnbCs4Lw+9jpmcnCeMo8Apmq609UV4wuBGFj9/DOEcDW0V+foEXOFbawG9XHePNnWe0lhJSHCypp8lkJTsligHxEVrL6TfC8HuYWTlJTM9ObJ/oFfiGU1VNHCqpxxrk+WzK6lt5Vynhk4NlWksJKZwBG8FyXwvD72Ey4iP4122TWTp9iNZSQoo3d53lxn/u4h9fndJaileZmpWIUS9xsKSOuhaz1nJChoJKewjtjFxh+AUCv2G7w7U2NStRYyXeJTrcwKTB8dhU2CHSN/iMZ26+iPXLZjErJ/AndkEYfq9gtansP1PHB3tLtJYSElQ0mMivbCLSqGP8oOCM6OmIs/7D1nxh+H3J4MRIIsP0WsvwCMLwe4HmNivfeVnht6uO09hq0VpO0OMc7U/JSgyJJHnOhVw7RTlGnxCMqdaD/y7RgNgIA+MHxWG1qewqEjent9nhMPwzQiSUc+yAWKLD9Zw+30xZXavWcoKeG/65i6uf3x402V5BGH6vMbM9fYN4HPcmqqqy3eHrnhEkERe9YdDruGXyIO6cNZQAX0Dq95TVtVJQ2cS5ulbS48K1luMxxAIuLzEzN5nnN50S2RS9zPmmNlrMNhKjjIxMj9Fajs948PIRWksICbZ1CBow6oNnnCwMv5eYkBlHVJiewqpmyutNQTVa8CdSYsLZfv9cSutbAz5VrsD/cEZOBZsbMXh+wvwMo17XXqFnl5iE8yo6ncSghMAset0fCquaeGPnGcrrhZ/fG6iq2n7vTg8yN6Iw/F5kalYiKTFhQRkV4A/YbCqNptCNmnpqfT7/8/kJNp8U7kRvcLamhdI6EwlRRoYHWd4tYfi9yG0zhrDlV3O4WR6ktZSg5GhZA9OWb+Le9w9pLUUTpjkWq4mwTu+w45S9X6cOTQg6N6Lw8XuRUIgp15LthdVYbCpRQbKoxl2c7ocdp0T9B29wxdh0UmLCiI0wai3F4wjD7wPqW8y0mG1igtfD7HSMyEI1FfOw1GiSo8OobGjj1PlmclKCyx2hNbERBhaMTNVahlfwyZBUluVEWZZflWX5sCzLtzi2/VKW5aWyLN/jCw1a8e9955j2p00892WB1lKCCrPVxh5HmcVgz8/THZIkMS3b/t13iDrPAjdw2fBLkjRJkqT/kyTpgCRJhZIkbZEk6TFJklxxYKcBPwQuA26RZXk2kKwoyutAoizL0/om3/8ZlhaNTRV+WE9ztLSB5jYrWclRIf0k9bWfXywU9CSfHy7n3vcPsSU/OCfOXXL1SJJ0OxAGPAc8AtQCSUAmcJ0kSQdUVf2qu/0VRTkBIMvyYOBZYDFwzPH2UcffO3vSUFFR4YrULqmtre3zvv0lRa8SHaajuLqFQwVnSY/9uhaslrp6IhB0fXHYfj1MyIjo17XhCbTsr2HxKpFGHaql7Rv9EAjn0Z/oqOvzA2dYl1fLiEQ9I+K0jcrzRn/1avglScoGPlNVtfPdVeH4t1eSpEGSJMWoqtrY3XFkWc4BlgOV2H84nEPgViCjNx1paWm9fcSr+/eHKVmlbMyrorBRz0W5F+rQUldP+LuuI5X2zKfzRg/0C61aaUhNVdn98KBuV5X6Q990hT/rUlWVg2UnALjkoiGkpWm/ItzT/dWrq0dV1VNdGP3OnynpyegDKIpSCFwCjAHMQJTjrVggOJ+nHDh90GIhl+f4/TWj+fMNY0MmP093SJIUVKkE/IHi6hbK600kRhkZFmTx+07cumIkSRosSdIVkiT9SpKkf0mSpLizv6IoNmAH8BYw3rF5DLDaneMEGs4JOGcUiqD/pMdF8K3xA0iKDuv9wyFAm8XGyYoex14CF3EO0KZmJQZd/L4TV1w9PwZ+gN1AhwOfAYeBT4Dfu9KILMvLHPtvBf6hKMpuWZYXyLJ8O1CrKMrmPuoPCEZnxBITrudMTYvI2yPwOOcb21jw9BbCDTp2PDgPfZAaK1+xs4PhD1Zcmdx9GLgFqMLuo48EXlZVtdjVRhRFeaaLbU+4un+go9dJ/OO7E8lOjiY5RoxQ+8tvPjmGBPxkbjYDEyK0lqM5yTFhpMSEUVLbyrGyBsYNjNNaUsCiqmr7k7nzST0YccXwX62q6mHH65skSboCWCVJ0qvAX1VVtXlNXRAhDw3ei8iXmK02Vh0spcVs4+cLcrSW4zdMz07iw33n2FFY7XeG32KxUFZWhslkQlVVbDYb9fX1Wsv6Bk5dz12ZitmqItWXke8HMl3pL0mSCA8PJyMjA4Ohd7Pe6yc6GH3n36slSfoS+A12182MXlsRCDzE4XP1tJht5KREkRorXGZOpmYn8uG+c+wuquXO2VqruZCysjKio6MZNGgQkiRhNpsxGv0vDUIg61JVldraWsrKysjMzOz1mL1O7kqSNFSSpIs6NWJSVfUR7L5/JElaKElSVJcHELTzh/+cYOHTW4KqhJuvcT6GB7P/tS84F3IpRTVYbarGai7EZDKRkJAgcgl5EUmSSEhIwGQyufR5V0b8RZIkzZQk6VngFFAKNOBYwCVJUgrwrKqqzf3QHRKcqWmhpLaV3UU1ZCaGXv54T+CMuAhm/2tfGBAfwZCkSIqrW/zOzx9ICeSqGk1EGvVEhekDRrMTSZJQVdd+9F0K51RVdRtwL1AAjAOWAMOAraqq/lJV1aI+ag0pnKOyXSKss0+YrTb2hnh+np5w9sn+M3UaKwlM2qwqZXUmiquD/4nc5eycqqpasYdwfuI9OcHNNEcWSZG3p28cr2ihxWwjNzWalBjh3+/MXbOzuGt2FkOSxNNkX2h2FEyKDg+80b67iCV/PmRkegzxkQZKaluFn78PJEYa+NHsodw4aaDWUvySoclRDE2OCnqj5S2azfYAxehw+3j4pZdeYtasWZjNZi1leQVh+H2ITie1h3WK9A3uk5kQzq8WDeeHs4ZqLcXvcdXXG4w88sgjvP7660yYMIHXX3+dBx54oNd9VFWluc1h+B2Ffe644w5yc3M5ceKEV/VqgbspG/4sSVJoJ0fpJ04/rEjfIPAGnx8u41t/28Er21xeXxl03HfffSxdupS0tDSWLl3Kww8/3Os+bRYbVpuKXi8R3qFyXn19PcePH/emXE1wtwLX28DlkiSlAseBdWooDy36wLzhydS1ZDN3eDLQprWcgCG/opF1B6u4fGK0qDTVA2aryonyRnaeruHq4f7p6x/3xKZu33t8yShuke1x6O8qZ/ntqm8a3ROPXdrj8ZOSLhybJiYm8thjj6GqKk899RSrVq2ipaWF3bt385vf/IY9e/aw7/Ax1qxdzyOPPYGUEQvA2rVrKS0tDUrD766r5xDwIWDCnlf/CUmSfiJJUqzHlQUp2SnR/NfCXCYNTtBaSkCx7lglz3x1jtd2nNFail8ztUM8v8XP4vm1JDs7mzlz5rBixQoARo0a1f7eiy++SHr6AEaNvYjzZWcBaGtr4+mnn+aFF14ISlePuyP+9dhz9bwNyKqqNkqSFIk9cdtCT4sTCJyI+H3X6BjPn1/VwsBeK134nsOPzHNphewtcmb76N8TGAwGDAYDNpvtgjmQhoYGJo4dwWULZlNWVgbAM888w1133cW4ceMoKgq+aHV3R/wNwAxVVZ92GP05gAU44HlpwUtlg4lXthXx3oFKraUEBG0WG3vPOOL3Rc6jXnGO+veVNGmsRDtOnjzJqVOnKCiw17rOz8/n4MGDyLLM888/z8qVKzl79iwWi4V77rmHG2+8kXvvvZfk5GRKSkrYu3cv1113HTqdDp1Oh9WqbRUuT+Nq6cXBwHexF1D5tSNcTAdcqqrqbOyLuwQucr6pjeVrTpIWY+RnlwbOqkatOFhST6vZRlZiuMhu6gJTsxL5YO859pWEbn7+4cOHk5+f3/73448/3v76gw8+uOCzF8tT2bZ9OzarFaPRSExMDO+88077+xs3bvS6Xl/jkuFXVfWMJEkbAD2wxbkZeNFbwoKZEWkxJEQZqWg0c6amhSFJIs1RT+xyFBKfNEj7EniBgHPEf7C0CYvVhkFU6OqR8vpW6lotDIgLI9kPk7R5A3dW7u4GdnfcJklSuscVhQA6ncTUoQmsPVbJzlM1wvD3wq7TdjfPpEEimscVBsRHcM/8bAZE2hDTuz2jqipNbVZQIUIfOk/ermTnXCs5fBGSJH0oSdIXjn9fAl95XWGQItI3uIaqqsRFGIgJ1zNxoBjxu8rPF+QyNyde1OPtBZPFhsWqYtBLhBlCp69cyc55WYc/nwc2OfL2IEnSZG8JC3Y6LuQKpOyFvkaSJJ69ZTxWm8r5KjEZLvAsTc78PI7VuqGCuz9xBmC4JEk5kiQ9A4i6d31keFo08RF6KhpMnD4vMlr3hqgj6x42m8qnR6t5+OMjWKyiSF53NJssAESFuxvZHti4a/i/BeQD72EP4ZzucUUhgiRJzM2J59JRqbRZhSe2O46VNmAyB1conS/Q6STe3FvBv/eVcrSsQWs5LrHqQCkLntrCqEfXs+CpLaw6UOrV9tr9+4gRf28cwx66maeq6iuASJPYDx5YkMnzt05gZLrwXXdFm8XGLS/uZsryTTS2WrSWE3BMdERBBUL9h1UHSnlk1THO1bWiAufqWnlk1bE+Gf8XXniBjz/+mLfeeqs9LPPxxx9n8+bNF3yuzfq1fz+8H/79mpoabr755i7fe/jhh3tNmHfgwAE++uijXtt59dVXefPNN/uksTPuftu1QBnwI4d/P/iWtAn8hoMldZgsNoYmRRITEVqP4p5g0kB7FFQgZIJ9akMBreYLXVKtZhtPbShw6zjvvvsuANdeey3f+c53OHz4MFu3biUrK+sbBjjcoGdURgxDkiL7NceWmJhIXNw3K56Vl5eTl5fHhg0betx/woQJXHfddb22k52d7bGFZO7eTcVAMuCc1BV3Yz9pbrOy70wtw9NiSBPFwy9A1NftHxMd4a9Kca3fx/OX1rW6tb07Vq1axYMPPtj+95QpU/j0008ZPXo069at4+mnn+amm24iPT2duro6du7cyR/+8Af+/ve/U11dzejRozl06BCqqhIVFcXKlSt57bXXqKqqoqysDJvNhtVqZePGjaxYsYJnn32WnJyc9tTGY3oAACAASURBVBXCHTl48CD/+Mc/ePjhh7n00kuprq7m9ttvZ9KkSURERFBQUEBERAQTJ04kIiICnU7HmjVrePXVV3niiSe44YYb+Oijj6iqqmLZsmXtx12/fn279ieffNKt/nHi7pWwAVgG3O74t6RPrQra+c0nR/nha/vYcFxErHTGOVKdKvLz9Im0mDCGJkXSZLL6vZ9/QHzXcSLdbe+Oznl4jEYjer3df79o0SLefvttXnjhBerr61mzZg133HEH69evp6GhgQkTJlBSUtKe0O2BBx5g2bJlbN++ncrKSq655href/99UlNTSUxM5KuvvsJoNHL11VeTk5NzgQ5VVcnLy2Pbtm2Ul5dz/vx5kpKSSExM5He/+x3f/va3GThwIM899xzDhg3DarVy6623Ul9fT1NTEzNmzCA9PZ158+ZhMBg4ffp0+7EbGhratfcVdw3/DlVVv62q6u2qqt4OfL/PLQsAmDJU5OfvijaLjX2O2rFTRH6ePjM1QOo833dJLhHGC81RhFHHfZfkunWcxYsXs3Xr1va/9+7dy5IlX49PIyIiGDhwIJOnTmP+VTfwg9vvwGw2I0kSl19+OVdffTVgT+gGcM0117By5cr2H4/y8nJmz57N3XffTVlZGdXV9lXlnV1FW7Zs4Qc/+AHf+ta3+MlPftLum3cet/NrJwsWLODxxx9n7ty5nDhxgvz8fEaPHn3Bj9mMGTNYunQpP/rRj/rs+nHX8EdLkvSUJEm/lSTpt8Bf+tSqoB3naHbX6ZqQrprUGad/f0R6DEnRIj9PX5k1LJlp2YlkuDly9jVLJgzgiSWjGRgfgQQMjI/giSWjWTJhgFvH+e53v4vVauXjjz/mgw8+YPTo0UybNo1hw4axadMmXn75ZR544AHefOtdzp49y5z5C7n88svZvHkzt9xyCyUlJe0J3QDCwsIYP34806ZNA+Dmm29m/vz5vPbaayxZsoQ9e/awYsUKysvLOXvWntK5pqaGF198kZYWe3nVuLg4nn/+efbs2cPJkycpLi4mPz+fEydO0NbWRn5+frur6Hvf+x4xMTEYjUYsFgurVq2ioKCAXbt2UVBQQF5eHu+++y5nzpxh0aJF7T9I7iK5Y2wkSVoGnAGqgVhgtKqqf+pTy27w0EMPqY899lif96+oqCAtLc2DijxDRUUFqampzPnfr6hsaOOzn01nWJr2ET7+0F8vby3iT2tPsnTaYB5ZPNJvdHWF0NUz+fn5DBs2rP1vs9nsUlpmb1Jc3Ux9i4WBCRHtAwt/0NUV7ujq3NePPvooy5cv/8bMtbuTs+OA4aqq/syRkrnczf0FnZAkiWlZiXx6qJydp2v8wvD7A7fPHMKi0aliRbPA46iqSpPJEb8fHlrx+07cdfUcBlY6Xm8FfuNZOaHJtCxH3h4/98P6EkmSGJwURWaif5YPDCQsVhsHz9Zx4GydJu1LkuRXbsxWR31dg14izI8jndzFndQv7n7rZiBbkqRZwAdAnpv7C7rA6ecvrzdprMQ/8CcjEQz850gFN72wmxUbCzVpPzw8nNraWr85r02ONA0x4YageaJUVZXa2lrCw10LCXfX1fMacCf2oixbgH+4ub+gC4YmRbLpvtl+PwHnK/626RSrj1Zwz/xsLh8jMn/3lylD7fWdlSJt4vkzMjLaI2BUVcVms6HTaTfSrm4y02qxYos00nr+a1eP1rq6wxVdkiQRHh5ORoZrtTbdNfxjVFVd4WhIAmbjQmpmWZZjgZexL/xarSjK3bIs3wlYgRTgfxVFCdlMUpIkCaPfgZ2na8grb0QXJKMxrcmIj2BoUiRF1S0cLW1gfGa8T9s3GAxkZn5dO1frSee9xbXszz/PzcMGXXDfaa2rO7yhy93Si1MlSdrv2KwDLsVu/HtjOnAb9qpd+2RZngLMVRTl+7Isfx+4CXjXTe1BSWOrJaTTE5jMVvY74vdlx0hV0H+mZSdRVF3CztM1Pjf8/sbFQxK4eEhoX1s+Kb2oKMo652tZlg8Di4GTjk1HgF/Qi+GvqKhwpakuqa2t7fO+3qSjroZWCz/5sIC6Vguf/HCMpqNdLftr/7lGTBYbuckRmBtrqehQNjYQzqM/0VHXqCS7q2DLiXKuGaFtxbdA6C9/whu63Cq9KEnSQWAS4FxRczPwlKvHcLh8irEXba93bG4FenVM9fdRxx8f4eBrXamqilktoK7VSq0axaj0WL/Q5Wvyjtgt/cxhqV1q8Pfz6G84dV0aEcfj685wsKyZpOQUzfP2aNVfK74sJDJMz3UTB3S5MNDfz6On6Euunnvpe66epcBvgUrAOeyIBc67eZygwx7PL8I6naUop4n8PB4lPS6CrOQowg06SmrdS3wWLFisNl7eXsSTa0/SZgnZKUWgb7l6bulLrh5Zlq8FPlYUpQF7euexjrfGAKvd1BGUOMM6Q7UOr8lsZf9Z4d/3Fm/cPpntD8xlaLK2rh6tOFLaQJPJytCkyJAPpnB3FjFakqSnAKfTaSxwS287ybJ8N3A/cF6W5TDgGWC3LMt3YHfzLHdTR1AyzZFQSymqwWZT0YVYuUFJknj6pnEUVjaTGCXy83ia1BBP++18kp6WnaSxEu1xNapH7yiwfgI4BDjLIbn0vKQoyt+Av/VJYQiRmRjJoIQISmpbOV7ewJgB3yzuEMyEGXRcOioNRmmtJLhpaLUQadRp7uf3NTtP2TNpCjei666eZZIkjQA+x56krdTxT5s14EHMtABJoysITH75wWGm/WkTB0vqe/9wENFmsbGn2O6oEIV9XHf1/By4Cujse8gGnvOoohDnu9MGc8W4dCaHWJxxo8nCo6uOMXtYMtdNFKWcvUV8pAGrTWVbYXVIxbIfOldPi9lGTkqUqHSH64b/KlVVj3TeKEnSBA/rCXnGDQwt946T3adr+PRQOSW1rcLwe5GZOUm8uess2wuruWd+Tu87BAlGvcSi0ankpERrLcUvcHUB1zeMvmP7Ac/KEYQq2wrt/tcZOWLizZtMzUpEJ8H+M3U0mSxEh4fGKvHxg+JZ8W0xTnUSWrM7AcLu0zX84t2DvLnrjNZSfMZ2h+GfKQy/V4mLNHLRoDgsNhWlyD9Xqgq8jzD8fkhlo4k1RytYdyw0CrBXNJg4WdFEpFHHhBDPI+MLnE9VzqesYKegsonNJ6toaetbfdpgRBh+P2R6dhKSBHuKa2k1B//F6hztT8lKJMwgLklvMzPEDP8He8/xozf2a1aPwB8Rd5kfkhQdxpiMWNostpB4HBduHt8yaXACf7lhHK98f5LWUnzC1gJ7RpjZw5I1VuI/CMPvp8zKtV+kzos2mJmQGc/07ERm5oob0xeEGXQsGZ9BSkzwhzVWNZo4Ud5IuEHHxYOFG9GJMPx+ysxc++h3a0HwP47fOiWTf902mZHpotC8wLNs6+BGDDeGZmH1rhCG30+ZPCSBCKOOE+WNVDaIWrwCz9JosvDIyqPc+tJuv6mF6w22OQZOs4Qb8QJCI4g3AAkz6LhtxhBiww0Y9MGbrO3zw+UMSohg7IDYkMsdoyVRRj3rj1dS02zm9PlmsoNwYZOqqmxxuEpnCf/+BQjD78fce8kwrSV4FYvVxm9XHaOh1cIX985iUEKk1pJCBp1OYkZOEp8fLmdbYXVQGv7qJjPRYQaIgRFpwff9+oMYYgk040hpAw2tFoYmRQqjrwHOKKrtQRrWmRwTxppfzOSze2YgaVjK1B8Rht/POVRSz3NfFnD6fLPWUjyO0/8q0jRog7Pfd5yqwWoLXj9/fKRRawl+hzD8fs7rO4tZsfEUX54IvlW8zogLZwSTwLdkJkYyJCmShlYLR84FV5rmNouN0rrQLDHpCsLw+zkz2+P5g+txvLnNyr4ztUgS7bWGBb4nWNM37DtTy/yntvCzt0Ueya4Qk7t+jtMPu7uohjaLLWhSGihFNZitKuMGxpEQJR7FteLKsekkRhqZE2RRL043YmaimDvqCmH4/Zy02HBGpMeQV97InuLaoPGH1zSbSYo2MntYcHyfQGVGTlLQXFMd2eKM3xduxC4Rhj8AmJ2bRF55I9sKqoPmJr1mwgCWXJSByeJS2WaBwGVqmts4UlqPUS8xZagos9gVweE3CHLa/fyFwZW3R6eTiAwTy+i1pqrRxKvbi3l9R7HWUjzC9sIaVBUuHpIgrq9uEIY/AJCHJJCZGMnItBhsQRB2d6a6mfoWs9YyBA7K6k38cXUeL20rCor0Dduc2ThF0r9uEa6eACAyTM+GZbO0luExfv+fPDbnn2fFt8ezcGSq1nJCnjEZsSRHh1FaZyK/sonhaYGdLG/X6RoAMX/UA2LEL/ApJrOVHaeqsdpUxg4IzcLy/oZOJ7Ubya9OBr47ceVPp/PP701kdEas1lL8FmH4A4i6FjOfHy7DYg3cCVGluJYWs41RGTGkxwV/PvhAYe7wFAA25we+4Y8M0zNveIpI09ADwvAHELe+pHDv+4fZfzZwV1ludowo5w5L0ViJoCOzcu3lPpWiGppMFq3l9JlgTj3hSYThDyBmO2KSN+VVaayk72x0aJ87XEy8+ROJUWFMGBSP2aoG7CremuY2Zj65mfs/PBwUk9TeRBj+AGL+CPsoeePJwDT8hVVNnD7fTHykgUmiDJ7fceXYNBaMTAnYpGZf5Z+ntsVMVWObcPP0gojqCSDkoYlEhenJK2/kXG0rAxMitJbkFgdL6pEkmDc8RRRd8UNumzmU22YO1VpGn3E+Cc8bIdyIvSEMfwARZtAxMyeJ9ccr2XSyilunZGotyS2unTCA2blJtJgDd3Ja4J9YrDa+ckxMzxeGv1fEsCvAcF7UXwRomuaUmHAGi8RZfkubxca2gvMBlwZ875k66losZCVHkZUcpbUcv8dnI35ZlucCjyqKcoksyzrgN0ABoFcU5V++0hHoLByZil53nFazDZtNRacLDF9mo8lCdJhe+F79HKWohttf28eItGgWBNDiuvXHKgC4ZFTgaNYSn434FUXZDDiHet8BShVFeQOYIcvyYF/pCHSSY8LYdv9cXr99csAYfYBfrzzKwqe3svNUYEaMhAry0ERiwvXkVTRxtqZFazkuoaoq64/bn1AWjRaG3xV87eNvc/y/GPi74/VJ4FLglZ52rKio6HOjtbW1fd7Xm/RHV0WjB4V0wtP9Zbba+OpkFU1tNsKtzVRU9C1OPBjPozfpq64pmTF8WVDHJ8opbpzgeX+5p/tLVVX+eOUQdhQ1MCDM1GdbEWznsSe0mtxNAWocr1uBjN52SEtL61eD/d3fW/RVl6qqHC9rJDc12ivFWTzZX5vy7EZ/RHoME4f1b0I62M6jt+mLrsUTrXxZUMf2s83cvcg738vT/ZWeDtNG9f84wXQee0Kryd1KwDkDEwsE/jpxH/PjN/dz7f/tbE9I5c+sOWofgV0+xj9vKsGFLBiRilEvoRTVUtVo0lqOwAtoZfg/B8Y7Xo8A1mukI2AZO9Ce4Gzdsb67wHyB2Wpjg8P/Kgx/YBAbYWBWbjI2FdYf8+/onjM1LVz/fzt5dVuR1lICCp8ZflmWLwJyZVkeB7wD5Miy/ENgq6Iohb7SESxc6ohe2HC80q/zk+w+XUNti5mclCiGpUZrLUfgIpeNSSMzMRK9nwcQrD9WwZHSBg6UBG7+Ki3wmY9fUZRDQMfonf/2VdvByJgBsWQmRnK2poU9xbVMzfLPEnPO3DyXjUkToZwBxDXjM7h+4gC/P2efHy4H7NeXwHXEAq4ARZIkFo9NB+CzQ2Uaq+meBy4bzr9uu5ibLh6ktRSBGxj0Or83+mdqWjhYUk9UmJ75w8VqXXcQhj+AWTzObvjXHqvw2xz9Br2O6dlJZIrVugFJeb2pfY7G3/iPY7S/cGSKqK3rJsLwBzCjMmLISo7CZoOiav9bbGMyW7WWIOgH9S1mFj69hf967yA1zW297+BjPj9sf9K90jEAEriOMPwBjCRJvPC9iWy5fw65fjZx2maxcckzW/nJm/sDurBHKBMXaWR6ThJmq8rqI/4VPXaqqoljZY3ERhhEUZ8+IAx/gDMkKQqjH6Y43lJwnsrGNkpqW4gSj+EBy5Lx9rWVqw761zzSgPgInr3lIu69JNcrCxiDHdFjQUJ9i9mvcqt8csBuKJaM9//IEEH3LBqVSqRRx57iWr+6viKMei4fk853p4o0X31BGP4gYMPxSmb/5Sv+tCZPaykANLZa2tNGX31Rr9k4BH5MdLihPePlp34cPSZwD2H4g4BxA+MwW218mVdFdZP2k3Brj1VgstiYmpUQcFXCBN/kW+MHALDyQKlf1LJdvjqP+z88TGFVk9ZSAhZh+IOA9LhwZg9LxmxV/WJU9onDH7zkogEaKxF4gpm5SaTFhjM0KYqmNm0jtVrNVj7Yd45PDpZhtmr/IxSoCMMfJFw/aSAAH+0v1VTH+cY2lKIawg06Lh8rVlMGA0a9jvX/NZP/++5EYsK1rda64XglDa0Wxg6MZWR6jKZaAhlRczdIWDgihfhIA0dLGzhe1sCojFhNdCTHhPHFvbM5VFJPfKRREw0CzxNu9I/IrH/vOwfADRMHaqwksBEj/iAh3Khvn0j9YO85TbWkxYaLEnhBiM2msq3gPNsLtamiVlrXytbCaox6iatE0EC/EIY/iLjxYvsoqKJBmxzq1U1t2Pw4U6igf6w5VsHtr+3jL+vyNWn/nd1nUVVYNDqNhCjxNNkfhKsniBgzII71/zWTwUlRvX/YC9z3wWGKq5v5683juWhQnCYaBN5jgcOdePhcPUfO1bfXhPAFFquNDxxunu9NE7H7/UWM+IMMrYz+6fPNbC+spqqxjaFJIiFbMBJh1HPNBHuk1jtKiU/bNuh1vH2HzP2LhnHx4Hifth2MCMMfhKiqyv4zdRwrbfBZm684KiBdfVEGcWJSN2i5dYq9ZvLHB0p9XpZxSFIUd87OEivBPYAw/EHIe3tKuOXF3fz1ywKftFfVaOLfjjDSO2YN9UmbAm3ISYnmklGptFlsvLHzrE/abDRZ/GLhWDAhDH8QcumoNIx6iY15VRRUen914xs7z9JmsbFwZIrfZQkVeJ47HT/ub+0+45PMqw/++wjf+tsODp6t83pboYIw/EFIckwY108aiKrC3zad8mpbTSYLb+0+A8Cds7O82pbAP7h4SAKXjkrlB9OH4O1x+LHSBtYfr6SouoWMeJH+w1MIwx+k/HhOFka9xGeHy7w66j9T00JMuIFJg+OZPCTBa+0I/Ivnb53Az+bneH0l7/ObCgH4tjyItNhwr7YVSgjDH6QMSoj0yah/VEYsa38xk2duushrbQj8G2+t3ThW2sC6Y5WEG3TiadLDCMMfxPhq1G/Q68RjeAhis6m8vqOYq57fTn2L2ePHF6N97yEMfxAzKCGSWyYP4vYZQ0mP8+yNU17fyvLVeX5Zi1XgGyQJ1hytoLCqmRe2FHn02LtO14jRvhcRhj/IeWTxSB68fLjHfbF//aKQV7YX88fV/lH8ReB7JEnigcuGA/CvHcWU1rV67NiNJgvJ0WH8eE6WGO17AWH4g5yOi12aTBZaPJBPfdfpGv69/xwGncTd83L6fTxB4DI+M57F49IxWWw89ulxj8XbLxyZyn9+PqM9dFTgWYThDxG2F1azeMV2Vmws7NdxGk0WHvroCKoKd83JIitZmxQRAv/hwcuGExth4Mu8Kj7c17/MsGarrf11fKTRb9JBBxvC8IcI0WF6yhtMvLStiM0nq/p8nOWr8yipbWXMgFh+OjfbgwoFgUpGfAS/WTwSgD+szutzUfY2i42lr+xhxZeFIsurlxGGP0QYnxnPPfNzUFX41YeHOdOHm/Pzw2W8v/ccYQYdT14/ljCDuHwEdr41PoPLRqcyNCmKtg6jdnf4/X9OsO9MHR/uO0eDD1YEhzLizg0h7p6bzYIRKdS1WPj5OwdobHXv5jp9vhmAX106jOFpouyd4GskSeIP147l3TunkJPiftqO9/eU8I5SQphBx3PfHi+qt3kZYfhDCJ1O4snrxzIkKZJjZY0sfXUPlW4Ubbl7Xg5v/nAy358u8qELvklshKH9KdBmU/n0UBkWF0b/r2wr4pFPjgHw6FUjGefDPP+hSkgUYlmbV8NLb+Rxrq4VCb6RXyQqTE+YXqKuxcKA+AjuuySXJY68405WHSjlqQ0FlNa1MiA+gqFJEew6XYtVtcczRxh0tJht6CWwqnzj/4HdHNfXxEUaeXnpJO58Yz/Hyxo4WtbAvG7C5aw2lRe3nGbWsOT2m1Eemtj+fuc+6fz9envf1WM88Z88ah0LhBIiDTxy5che2wF4akMB5+pa289BQqQBkKhtMXd5XjoeJz7SiMliocX89dWSEGnkkStH9HgOO+vtiPPa69w20Ot39ATd9bUrfdzbMTrzP5+f4K3dZ3ktM44nrx/XZRCAxWpj+ZqTvL7TnuvpgcuGc+PFgzz6nV1l1YFS/rw2j4rGgwyIj2De8GQ2nTzf6/fsqU+/vpYMNLdZabParyUJ+6K03y0Z3f65jtepL2yFpGW6U1mWfwlUAPGKoqzo7nMPPfSQ+thjj/WpjVUHSvn1J0cxWVz/nhFGHU8sGd3e8asOlPLIqmO0mvvmu+zuuBUVFaSlpfXrmH2luqkNpaiWy8bY27fZVLYWVpMRF05tTTVHq2HlwVKOnGsgNzWaT346DYP+6wfErvqk4/fr7X1Xj/HwyqOYrReeO4MOll87ttt2jHoJVVWxuHi6Iow6rpswgI8OlPZ6jo16iT9eM+aCm9J5HrvT29vxrDaVznOZHb9jX+l4fXXX19dNGMAH+8712MdOXDmnTpSiWn75wSHK6k1EGnXcMGkgC0amkpUchbmpluzMAaiqyq0vKRw+V88frx3LkvHa1NF15f7u6nv21KeuXEszshPYd7a+y891bK8/duLRRx9l+fLl3yhgoJmrR5bl2UCyoiivA4myLE/zRjtPbShwy+gDtJptPLXh61z2T20o6LfR7+q4WpIUHdZu9AHWH6/kztf3cfXzO/jeW3n8YXUeR841kBYbzn9fMeICow9d90nH79fb+64eoysjarHRYztmq+tG39nme3tKXDrHZqva7TnsTm9vx+sqgKXjd/QE3fX1e3tKeu3j3o7RlU55aAKr7p7OkvEZtJhtvLHrLHe8vo9LntnK9tP1gH1e4MHLh/Ov2yZrZvTBtfu7q+/ZU5+6ci1tP1Xb7ee8bSu0dPUsBo45Xh91/L2zuw9XVFT0qZG+riYsrWttb9OTKxI7Hre2ttZjx+0vp8vOc/GgaKqaLDSaLEwcFMOsrDhmZccRHWb9Rv931yfO79fb+/05hqufcQd37HXH7wBfn0dPXiddteMuHa+v7rT19L07t+/KOe3MA3PSWDIihi2n6tleVE9ti5WGpub2zw9yeBn78z37i6vnzdX+cPO3v9f2vGEntDT8KUCN43Ur0ONPfl8fdQbER3CuDzfkgPiI9jb7eozejgt9/16e5q6Fady10P7alUfL7vrE+f16e78/x3D1M+7g9K+6QudzCHhUS0/tuEtvfd3T9+7cvivntGsNMK9D8lYtXZxd4ep5c7U/3LmWXG3P0/2lZVRPJeCc8YkFznujkfsuySXc4F6Nzgijrn3SzXmMCGP/u6rzcQOZrvqk4/fr7X1Xj2HUf/PcGXT02I5RL+HOEoMIo46bJw9y6Rwb9VK357A7vb0dT9fFLh2/oyforq9vnjyo1z7u7RiBfk27cn939T176lNXrqUZ2Qndfs7b/arliP9z4ErgPWAMsNobjSyZMIC6hnpe2lXZ56ge5+tgiOrxFF31Scfv19v77hyjp4iT7o7h3OZOVM/kIQn9iurpSm9HtIzq6amvJw9JcKl9V85pIOLUb4/qMbsc1dNbn4qonm6QZfkRoARIVBTlqe4+15+oHvC/R0snQpd7CF3uIXS5RzDq6i6qR9M4fkVRntCyfYFAIAhFxMpdgUAgCDGE4RcIBIIQQxh+gUAgCDGE4RcIBIIQQxh+gUAgCDGE4RcIBIIQQxh+gUAgCDECJh//o48+qrUEgUAgCAo0XbkrEAgEAt8jXD0CgUAQYgjDLxAIBCGGMPwCgUAQYgjDLxAIBCGGMPwCgUAQYgjDLxAIBCFGwMTx94Ysy78EKoB4RVFWdNg+ArgFaAZWKYqS19U2H2u6FVgGxAFLFUVRHNufA24C9iuKcoU3NPWky/Hex8B07P3yI1/1VU+6ZFleA4zCXsBKVRQluyutXtQ1F3hUUZRLOm2fCczCPoB6RVGUiq62aaBrGfB97P11o6Iopxzbte4vHbADGAL8TVGUx7XuL1mWJeAQEOPYVKwoytyutHpJUyzwMjAZWK0oyt0d3ksD7gHKsNuEbV1t60u7QTHil2V5NpCsKMrrQKIsy9M6vP1X4GlgBbC8h20+0eS40JoVRZkG/AV4zLF9ELBXUZQMLxv9bvtKluUpwN8dGpyGwet91ZMux41xn6IoQ7Eb/7d70OoVFEXZDER28dYfsZ/Dt3Gcx262+UyXLMsJ2A3CxdjLm/7Ksd0f+ut64PsODU5Dqml/AYOByxRFyQJmAx/3oNUbTAduA8YBlzjOk5M/AG8oivI34GGH7ehqm9sEy4h/MXDM8fqo4++dsixHArmKojQCyLKc7TAknbcZFEWx+EKToigqsNKxfTcww/F6IfAbWZZvAG5TFKXKw3p61OX4ewHwc1mWvwB+in3E6Iu+6laXoigNwBHH9suAtV1pVRSl2QuaOtLW8Q/Hk5DFcT6LZVme09U2L2v6hi5FUWqBjY4/dwPjHa817S8Hs4DnZFl+E3gAGIb2/VXc4c/rgI+60qoois0bYhRFWed8LcvyYewjeSeXAR1/pLO62XbK3XaDYsQPpAA1jtetQIbjdSJQ3+FzFuzulc7bUn2oqSOXAk8BOEa6ucAG5zYv0a0uRVGeBLKBKuAhuu4/b/RVj7o6MAf4qhutvqajXrD3+IH0fwAAA8JJREFUVVfbtGQG8E/wi/5CUZR7sV/jmdhHuf7WX9lOt1gXWr2KY0BarCjKmQ6bjY4fRfj6nuhqm9sEi+GvBKIcr2OB847X54GIDp+LAhq72FbrQ00AyLI8DChSFOWoc5uiKKqiKE8DYV7Q45Iux2j+QexGoqv+80Zf9apLlmUDYFUUxdqNVl/TUS+AqZttmuBwGazt6DPXuL+cGpqBXwAT8a/+SgXKO27rpNXbLAV+22lbY4fXznuiq21uEyyG/3O+fqQdA6yRZTleURQTUCTLcpQsyxHAGUVR6rrY1uIrTQCyLKcDExRF+VCW5RhZlqOdvjpZlsOwP6J7i550Of2FscCWbvrPG33Voy4HC4AvnX901uolTd9AlmW9LMuxiqKcxPGjKMtyDrCxq22+1uV4PRxIUBTlS1mW02RZlrTuL8drp4ZkYIO/9JeDa/jaBfsNrV7Wci3wsaIoDbIsp8uy7Hzy2egYIAKEOwIrutrmNkGTpE2W5UeAEuyPixuAhxRFuVWW5XHAjdhHEysVRTna1TZfacI+I78Bu9sEQAJk4D2gCbu//TVFUZq8oak7XY6+2grsc/x7VVEUq6/6qiddjvf+gD0iw+z4+xtavajrIuw/TFdiHy3PVhTlQVmWF2I/dxHAC4qilHa1zZe6sLt2PgcasF9bpYqiXK11f2GfuN0BrAb2KYrinKTXtL8URXnQ8d5fFEVxToRHdaXVS5ruBu7HPnIPA94ArnZEFg0Efo7d779PUZTNXW3rS7tBY/gFAoFA4BrB4uoRCAQCgYsIwy8QCAQhhjD8AoFAEGIIwy8QCAQhhjD8AoFAEGIIwy8QuIEkSeskSRL3jSCgCZZcPQKBR5Ek6Vns6QRGAmnY1xUsAC5XVdUreVsEAl8h4vgFgi6QJGm0qqrHJEm6DRilqupDzm1aaxMI+osY8QsEXdCNgTdKkrQFmA/8CXumx0hgLPA89roFR1RVfUKSpKnYs08uAlaqqvpxF8cTCDRB+CoFAhdRVfWg438L9hTSpaqqLsOeuK5QVdVbgcsdH/8VUA18gf2HQSDwG8SIXyBwD0uH/50pq5s6vHbmwBmjqupqADEZLPA3xAUpEHgHVZKk6x2vF2uqRCDohBjxCwTdIElSAjATyJUkKRN79sRcSZJygUnYjXsmMByYLklSJJAjSVIO9iysr0mSdC/wA22+gUDQNSKqRyAQCEIM4eoRCASCEEMYfoFAIAgxhOEXCASCEEMYfoFAIAgxhOEXCASCEEMYfoFAIAgxhOEXCASCEOP/AUrgE+BxX4MaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1)\n",
    "ax.plot(times, true_intensity_function(times), \"--\", label=r\"True $\\lambda$\")\n",
    "ax.set_xlabel(\"Time\")\n",
    "ax.set_ylabel(\"Intensity ($\\lambda$)\")\n",
    "ax.scatter(arrival_times, torch.zeros_like(arrival_times), label=r\"Observed Arrivals\")\n",
    "ax.legend(loc=\"best\")\n",
    "None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameterizing the intensity function using a GP\n",
    "\n",
    "When using a GP to parameterize an intensity function $\\lambda(t)$, we have to keep two concerns in mind:\n",
    "\n",
    "1. The intensity function is non-negative, whereas the output of the GP can be negative\n",
    "2. The intensity function can take on large values (e.g. 900), whereas GP inference works best when latent function values are fairly normal (e.g. $\\mathbb{E}[f] = 0$ and $\\mathbb{E}[f^2] = 1$).\n",
    "\n",
    "As a result, we will use the following transform to convert GP latent function samples into intensity function samples:\n",
    "\n",
    "```python\n",
    "# function_samples = <samples from the GP>\n",
    "# self.mean_intensity = E[ \\lambda (t) ] - computed from data\n",
    "intensity_samples = function_samples.exp() * self.mean_intensity\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Cox process likelihood\n",
    "\n",
    "Given an intensity function $\\lambda$, the likelihood of a Cox process is\n",
    "\n",
    "$$\n",
    "\\log p(t_1, \\ldots, t_n | \\lambda) = \\left[ \\sum_{i=1}^n \\log \\lambda( t_i ) \\right] - \\int_0^T \\lambda(\\tau) d\\tau\n",
    "$$\n",
    "\n",
    "The first term, which are the arrival log intensities, is easy to compute, given the `arrival_times`:\n",
    "\n",
    "```python\n",
    "# arrival_intensity_samples = <a function involving arrival_times>\n",
    "arrival_log_intensities = arrival_intensity_samples.log().sum(dim=-1)\n",
    "```\n",
    "\n",
    "The second term, which is the estimated number of arrivals, can be computed analytically.\n",
    "However, it is much easier to compute this integral using quadrature.\n",
    "Using a grid of `quadrature_times` evenly spaced from $0$ to $T$, we can approximate this second term:\n",
    "\n",
    "```python\n",
    "# quadrature_intensity_samples = <a function involving quadrature_times>\n",
    "# max_time = T\n",
    "est_num_arrivals = quadrature_intensity_samples.mean(dim=-1).mul(self.max_time)\n",
    "```\n",
    "\n",
    "Putting this together, we have\n",
    "\n",
    "```python\n",
    "log_likelihood = arrival_log_intensities - est_num_arrivals\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Using the low-level Pyro/GPyTorch interface\n",
    "\n",
    "Unfortunately, this is not a likelihood function that is easily written as a GPyTorch Likelihood object.\n",
    "Because of this, we will be using the low-level interface for the Pyro/GPyTorch integration.\n",
    "This is the more logical choice anyways: we're simply trying to perform inference on the intensity function rather than constructing a predictive model.\n",
    "\n",
    "Here's how it will work. We'll use a `gpytorch.models.ApproximateGP` object to model the non-homogeneous intensity function. This object needs to define 3 functions:\n",
    "\n",
    "- `forward(times)` - which computes the prior GP mean and covariance at the supplied times.\n",
    "- `guide(arrival_times, quadrature_times)` - which defines the approximate GP posterior at both arrival times and quadrature times.\n",
    "- `model(arrival_times, quadrature_times)` - which does the following 3 things\n",
    "\n",
    "  - Computes the GP prior at arrival time and quadrature times\n",
    "  - Converts GP function samples into intensity function samples, using the transformation defined above.\n",
    "  - Computes the likelihood of the arrivals (observations). We will use `pyro.factor` rather than a GPyTorch likelihood or a `pyro.sample` call (given that we have defined a simple function for directly computing the log likelihood).\n",
    " \n",
    "Putting it all together, we have:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "class GPModel(gpytorch.models.ApproximateGP):\n",
    "    def __init__(self, num_arrivals, max_time, num_inducing=32, name_prefix=\"cox_gp_model\"):\n",
    "        self.name_prefix = name_prefix\n",
    "        self.max_time = max_time\n",
    "        self.mean_intensity = (num_arrivals / max_time)\n",
    "        \n",
    "        # Define the variational distribution and strategy of the GP\n",
    "        # We will initialize the inducing points to lie on a grid from 0 to T\n",
    "        inducing_points = torch.linspace(0, max_time, num_inducing).unsqueeze(-1)\n",
    "        variational_distribution = gpytorch.variational.CholeskyVariationalDistribution(num_inducing_points=num_inducing)\n",
    "        variational_strategy = gpytorch.variational.VariationalStrategy(self, inducing_points, variational_distribution)\n",
    "\n",
    "        # Define model\n",
    "        super().__init__(variational_strategy=variational_strategy)\n",
    "\n",
    "        # Define mean and kernel\n",
    "        self.mean_module = gpytorch.means.ZeroMean()\n",
    "        self.covar_module = gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())\n",
    "    \n",
    "    def forward(self, times):\n",
    "        mean = self.mean_module(times)\n",
    "        covar = self.covar_module(times)\n",
    "        return gpytorch.distributions.MultivariateNormal(mean, covar)\n",
    "\n",
    "    def guide(self, arrival_times, quadrature_times):\n",
    "        # Draw samples from q(f) at arrival_times\n",
    "        # Also draw samples from q(f) at evenly-spaced points (quadrature_times)\n",
    "        with pyro.plate(self.name_prefix + \".times_plate\", dim=-1):\n",
    "            pyro.sample(\n",
    "                self.name_prefix + \".function_samples\",\n",
    "                self.pyro_guide(torch.cat([arrival_times, quadrature_times], -1))\n",
    "            )\n",
    "\n",
    "    def model(self, arrival_times, quadrature_times):\n",
    "        pyro.module(self.name_prefix + \".gp\", self)\n",
    "        \n",
    "        # Draw samples from p(f) at arrival times\n",
    "        # Also draw samples from p(f) at evenly-spaced points (quadrature_times)\n",
    "        with pyro.plate(self.name_prefix + \".times_plate\", dim=-1):\n",
    "            function_samples = pyro.sample(\n",
    "                self.name_prefix + \".function_samples\",\n",
    "                self.pyro_model(torch.cat([arrival_times, quadrature_times], -1))\n",
    "            )\n",
    "        \n",
    "        ####\n",
    "        # Convert function samples into intensity samples, using the function above\n",
    "        ####\n",
    "        intensity_samples = function_samples.exp() * self.mean_intensity\n",
    "        \n",
    "        # Divide the intensity samples into arrival_intensity_samples and quadrature_intensity_samples\n",
    "        arrival_intensity_samples, quadrature_intensity_samples = intensity_samples.split([\n",
    "            arrival_times.size(-1), quadrature_times.size(-1)\n",
    "        ], dim=-1)\n",
    "\n",
    "        ####\n",
    "        # Compute the log_likelihood, using the method described above\n",
    "        ####\n",
    "        arrival_log_intensities = arrival_intensity_samples.log().sum(dim=-1)\n",
    "        est_num_arrivals = quadrature_intensity_samples.mean(dim=-1).mul(self.max_time)\n",
    "        log_likelihood = arrival_log_intensities - est_num_arrivals\n",
    "        pyro.factor(self.name_prefix + \".log_likelihood\", log_likelihood)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = GPModel(arrival_times.numel(), max_time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Performing inference\n",
    "\n",
    "We'll use Pyro to perform inference. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Defining the quadrature times\n",
    "\n",
    "We'll use 128 quadrature points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "quadrature_times = torch.linspace(0, max_time, 64)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Pyro SVI inference loop\n",
    "\n",
    "Now we'll use Pyro to perform inference."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8c7c4baf232f4820b0133825fc729814",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(IntProgress(value=0, max=200), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# this is for running the notebook in our testing framework\n",
    "import os\n",
    "smoke_test = ('CI' in os.environ)\n",
    "num_iter = 2 if smoke_test else 200\n",
    "num_particles = 1 if smoke_test else 32\n",
    "\n",
    "\n",
    "def train(lr=0.01):\n",
    "    optimizer = pyro.optim.Adam({\"lr\": lr})\n",
    "    loss = pyro.infer.Trace_ELBO(num_particles=num_particles, vectorize_particles=True, retain_graph=True)\n",
    "    infer = pyro.infer.SVI(model.model, model.guide, optimizer, loss=loss)\n",
    "\n",
    "    model.train()\n",
    "    loader = tqdm.notebook.tqdm(range(num_iter))\n",
    "    for i in loader:\n",
    "        loss = infer.step(arrival_times, quadrature_times)\n",
    "        loader.set_postfix(loss=loss)\n",
    "        \n",
    "train()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluate the inferred intensity function\n",
    "\n",
    "Here we'll see how much our inferred intensity function matches the true intensity function. We'll plot the intensity function that corresponds to the GP mean, plus the intensity function at the 5th and 95th percentiles of the GP."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here's a quick helper function for getting smoothed percentile values from samples\n",
    "def percentiles_from_samples(samples, percentiles=[0.05, 0.5, 0.95]):\n",
    "    num_samples = samples.size(0)\n",
    "    samples = samples.sort(dim=0)[0]\n",
    "    \n",
    "    # Get samples corresponding to percentile\n",
    "    percentile_samples = [samples[int(num_samples * percentile)] for percentile in percentiles]\n",
    "    \n",
    "    # Smooth the samples\n",
    "    kernel = torch.full((1, 1, 5), fill_value=0.2)\n",
    "    percentiles_samples = [\n",
    "        torch.nn.functional.conv1d(percentile_sample.view(1, 1, -1), kernel, padding=2).view(-1)\n",
    "        for percentile_sample in percentile_samples\n",
    "    ]\n",
    "    \n",
    "    return percentile_samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAEFCAYAAAABjYvXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOydd3ib5dX/P7e2LMt7xc7em0AUCHuvpNBCKaUUWnih0AW8FNpfJzRtaUNpeQsECoWOl+4GKJTS8jJC2SNKyCTTiffWsGRt6Xl+f0h2LFteiTUS35/r8hVbesaxI+nc5z7nfI9QVRWJRCKRTGw02TZAIpFIJNlHOgOJRCKRSGcgkUgkEukMJBKJRIJ0BhKJRCIBdNk24HD55je/KcugJBKJ5DBYu3atGPjYUesMANasWXPY53Z0dFBRUTGO1owP0q6xkYt25aJNIO0aK8eqXXfffXfKx+U2kUQikUgyFxnYbLYzgLvtdvu5NptNA3wPqAW0drv9f1M9linbJBKJZKKTscjAbre/AZgTP14NtNrt9j8AJ9tstilDPCaRSCSSDJDpnEE48e8q4JeJ7/cB5wHnp3jst8NdrKOj47ANcbvdh31uOpF2jY1ctCsXbQJp11CEQiFCoRCKoiQ9rihK1m1LxVjsMhgM5OXljerYbCWQywBX4vsgUDXEY8NypMmdXEwOgbRrrOSiXbloE0i7UlFbW8vkyZMxGo0IcajIJhKJoNfrs2bXUIzWrlgsRl1d3aj/ttlyBp1Ar7uyAo4hHpNIJJK0oqoqJpMp22aMO1qtlrEIkWarmuhfwNLE93OBV4Z4TCKRSCQZIGPOwGazLQFm2Wy2xcBfgJk2m+2/gLftdvuBIR6TSCQSSQbI2DaR3W7fDvSvEPr2gOdjAx+TSCQSSTKBiILQKOi047uWl01nEolEchQRiCjEVJVf//rXnHrqqUQikXG5rnQGEolEkkG++93v8vvf/57jjjuO3//+93zjG98Y9bmqqhJT4knhG264gVmzZrFnz55xseuo1iaSSCSS8WLe3eNXs7JnzXlDPve1r32NkpISnnzySa699lpcLteQxw5EUUHpVyHk8XjYvXs3ixcvPiJ7QToDiUQiySglJSVJPxcXF7NmzRpUVeX+++/n+eefJxAIsHHjRr73ve+xadMm9u3bx6uvvsoP7/kxijZegf/SSy/R2trK7t27x8Uu6QwkEomEwav5TDadzZgxg5qaGmbOnAnA/Pnz2bhxIwBPPPEEV111FcuWLePgwTqqZy8kHA7zP//zPzz++OPcd99942KDdAYSiUSSA+h0OnQ6HYqiJDWLeb1eZs6cyemnn86+A/WoqspDDz7ATTfdxOLFi6mvrx+X+8sEskQikWSYffv2cfDgQWprawHYv38/27Ztw2az8fDDD/Pcc8/R1NRENBrlq1/9KldccQW33nor1sIi2lpb+XDzh1x22WVoNBo0Gg2xWOyIbZKRgUQikWSYOXPmsH///r6ff/CDH/R9/9RTTyUdu3LlSt5//30AunpCVE2axB/+9Ke+5//zn/+Mi00yMpBIJJKjhKiiJonpjSfSGUgkEslRgKqqRGMq6XEF0hlIJBLJUUGi14w0BQbSGUgkEsnRQG/ncbqQzkAikUiOAqQzkEgkEol0BhKJRCKBqKKkLV8A0hlIJBLJUUEkpqJJ4QzGS8paOgOJRCLJIOvXr+fMM8/kwQcfHLV89csvv8wD9/8MkaKwdLykrKUzkEgkkgQmk6nvy2q1Jv38xBNP9B33xBNPJD3X/2skVqxYwZw5c7j11lvZu3cvW7duHfGcOXPmEAqHh9wm6pWyPhKkHIVEIpFkiWAwyLp166ipqSEvL49zzz23T6567dq1uN1uNmzYQJfDMeQ1xkvKekI6g/7TgiQSiaSXYDDY9/1wEtY33ngjN95442Hfp6Ojg6eeeoqbb74Zr9dLTU0N5557Ll/60pf65Krr6up47LHHuP/++2lua+eJ3z456DrjKWU9IbeJGpwBXtnrzrYZEolkglJRUcEVV1zBZZddBsTlq+GQXPWXvvQlysvLcblcOJ1O1CHWrr/4xS/GTcp6QjqDqKJS6wiiyOhAIpFkmA8++ICDBw/iSGz99MpXA0ly1SUlJdx6663cfvvt/Pvf/6a1pZnXN7zKP599BoDm5mY2b948blLWE3KbCMAdjODwhSm3GrNtikQimUBceeWVXHnllX0/95ev7i9XDXD66adz+umn4/SF+XQkhlYjCEfjH/g1NTX85S9/6Tv2SKWsJ2RkAOALKbR6giMfKJFIJFkmndLVvUxYZxBWVA50+LJthkQikYxINJbe7mOYwM7ArNewt6Mn22ZIJJIsI4QgGAwmzR3OJRRVRVEZ8xyDWCw2pmhiwuYMzHoNLn+EnmCUfNOE/TNIJBOe8vJy2tvbiUajSQ5BURQ0muyvlxVVxROIok1oUUQVBb9J3/fzUAghKCwsHPV9JvinoKDVE2SOKT/bhkgkkixRUFBAQUHBoMc7OjqoqKjIgkXJ1Hb6eHpnHdWF8e7mg+1ubjt/OtNK88b1Ptl3e1lEAPUOf7bNkEgkkiHxBCOQgR2sCe0M8k1amTeQSCQ5jaMnnFKtdLyZ0M7AYtDS5AoQjirZNkUikUhS0u4NYdSl/6N6QjsDTcLdtntDWbZEIpFIUtPpDWHUa9N+n6wlkG02Wx7wLWAzcBLwY+BKIAaUAT+32+1pX7IrCjS7AkwpNqf7VhKJRDImVFXF4QtTajGk/V7ZjAwuBLrsdvvfgUbgBuAMu93+W6Ad+FQmjMgzaNgnm88kEkkO4g/HiMbUEctIx4NslpZ+AKyx2WwvAFbAD+xLPLcTuBX463AX6OjoOKwbO5xBwqEwPp8PEVPZ0dhJ2wwDmnS3+I0Ctzs31VSlXaMnF20CaddYyQW7OnrChIJBfL5DAnThcBiHw4E5Nr7FL1lzBna7vdlmsz0A/Ar4PTAJcCWeDgJVI13jcGuAvcKHwdiJxWIBwN8dRGcpoiw/N0TrcqG2ORXSrtGTizaBtGusZNsul9qDyezBYjn02WToiVBaWkrFsdJnYLPZpgA1wMXA5wE90PvbWYGhR/uMNyq0dEvROolEklt4AmFiamaqHbOZM1gOuOx2ewj4BfG2ikWJ5xYCL2bKEJ1WSNE6iUSSc7R7whi0mfmYzqYzeBGYYrPZVgHzgIeAjTab7QbiW0Z/zJQhVqNONp9JJJKco80TxJSBHgPIbs4gCHwj8eO/Ev+uy4YtJr2G1u6QFK2TSCQ5Rbs3hCkDPQYwwZvOeonLvAo57EYikeQM4aiCNxhFr81MlaN0BgmEkKJ1Eokkd3AHImgEaZ9w1ot0BgmsRi272rzZNkMikUgA6A5EMno/6QwS5Bm1tLiDBCOxkQ+WSCSSNOP2h4llcPqadAYJNEIgELLfQCKR5ARtnlDGykpBOoMkVKHKvIFEIskJWrszV1YK0hkkYTXqZN5AIpHkBB3ecMbKSkE6gyTiw26ChGTeQCKRZJFwVMEXylxZKUhnkIRGI1BVlVaPHHYjkUiyh8sfQQiRsbJSkM4gJTJvIJFIskm8rDRzlUQgncEg8o1adsu8gUQiySIufxgls75AOoOBWIw6GlwBwtHMyMZKJBLJQNq6QxgymC8A6QwGoU3kDdqkTpFEIskSrZ4gRn1mP56lM0iBokKDU+YNJBJJduj0hjDpMldWCtIZpCSeN5DzDSQSSeYJRmL0hGMZLSsF6QxSkm/UcdDhJxKTeQOJRJJZugMRtBkuKwXpDFKi1QgURaVd9htIJJIM4w5Es3Jf6QyGQAUaXYFsmyGRSCYYLl8YJYNqpb1IZzAEeQbZbyCRSDJPmyeY8bJSkM5gSKxGLQe6/ERl3kAikWSQ1u7MzT3uj3QGQ6DTaojGFDq84WybIpFIJhCd3hCmDPcYgHQGI9Lokv0GEokkMwTCMfyRGDqN3CbKKcx6DXvaZb+BRCLJDN2BCBpBxstKQTqDYbGadOzv9BHLtGKURCKZkLgDkUyLlfYhncEw6LQaIlGFTq/sN5BIJOnH5QuTrZIV6QxGQZPsN5BIJBmgpTuEUZf5LSIAXVbuehRh0mvY1d6DbXpxtk2R5ABOX5htTd28e9BFeb6BE6cXM7sinzxD5ksBJcce7d5gxgXqepHOYASsJh37O3pQFBVNFjL8kuwTjirs7+zh7VontZ0+hIBis56W7iB/2tiERsCi6gJs04qZW2HJtrmSoxRVVenwhigwZedjWTqDEdBrNURiCp09ISoLTNk2R5JhVFXlr5ua2N7sId+oZVKhEU2i0sNs0FKcpyemqOxt72FLUzefWDqJOQVZNlpyVBKIKIQiCnpLdnbvZc5gFCiqzBtMVHa2eNje5GFykYniPEOfI+iPViMoyzdQZTXy3NZWdnfI3hTJ2HH7I2ShorQP6QxGgUmvYbfsN5hw9ISiPLOlldJ8w6jqvg06DWX5Bv6+3UG9QzoEydhwB8JkQZ+uD+kMRoHVqGN/hw9F9htMGFRV5d872glGYmNKDpsNWiwGDb99t0GWJEvGRIcnhCB7oUFWcwY2m00Anwc6gK3AVYnvC+12+7ps2tYfg05DyBemyxemwmrMtjmSDFDb6WNjvYuawrHnifKNWsLAb95p4MtnzsCapYSg5Oii0RXAZMje+jzbkcFaYKPdbv8XMAMotdvtvweKbTbbSdk1LRmZN5g4BCMx1m9upsisP+wKsmKLHm8wwh/eb5QRpWRUNLkD5GVBoK6XrC1ZbDbbycBJQJPNZruGeBP2rsTTHwGrgPeHu0ZHR8dh3dvhDBIOhfH5fKM+R41E2bS/lcmm9KqYut3utF7/cJlIdr2y10WHq4dKqx6fb+xbPcFgEACLRmVXs5+3PxLMq8gbbzPHzET6PxwPMmlXKKrQ7uqhMl9HNDz8AiQcDuNwODDHxjePmc349RPAb+x2+5M2m+0x4Hrg8sRzQaBqpAtUVFQc1o29wofB2InFMvqacL1RodUfo7y8PO0iUof7e6WbiWBXoyvA1o5OplUUoj2CvpLe11aVLsq7LRFOXlCGTpvtQHxi/B+OJ5myq8UdJM/sJD9/5G1oQ0+E0tJSKkrHd4Ex6lenEOJ4IcSjQoitQogDQoi3hBBrhBA1h3lvE+BJfP9P4I9A729nBRyHed20YNBpCERjdPXI+QbHMht2d2LSa47IEfTHatLh6Amxrdkz8sGSCYvDF0bNZikRo3QGQojrgROBh4BzgbnEV/HPAZcJIU4/jHu/BRyf+F4P7AeWJn5eCLx4GNdMLyq0dAezbYUkTbR7Quxq81Ji0Y/rdYvz9Px7ZzvhqJyaJ0lNmydItgUORnQGQogZwAuqqj6mqupOVVW7VFWNqqraoarqZlVV1wEHhBD5Y7mx3W5fD1hsNtungWnAfUDQZrNdD7jtdvsbh/H7pBWDTsMeORf5mOWdWgdajUjZWHYkWIw6PIEImxtc43pdybFDg9OPOcv6ViPmDFRVPTiKY5oP5+Z2u/3OAQ/96HCukykKTDr2tPtQVTUrwyck6aM7EGFjvSttpcOl+Qb+76MOlk0pysp8W0nuoqpqPGeQZWcwpoyWEGKKEOIiIcSdQoj/FULY02VYLmLQafCHozh8Mm9wrLGxzoUK45YrGIhZryUQUXjvoIwOJMn4wzH84Rh6bXYXmKPZJrpZCPGOEMIN7AVuBPKBfwBXp9m+tBCMxA47WaMCDU7Zb3AsEYzEeHO/gzKLIa33Kc83sGF3Jz2haFrvIzm6cPrCCCGyvtswmsjgW8DtwHLiVT8m4Deqqj6tquredBqXDjY1uPnqX7ZR5zo8qQCzXsPOVpk3OJbY0ugmFFUw6NJb+mnQxRVw396fU4Vykizj8EVQslxJBKNzBh9TVfV9VVVrVVX9FLAOeF4IcbsQIvuF02PE7Y/Q0h1kT0eA7kBkzOcXmHTsbe8hGpOVIccC0ZjCq3u6KM4b3wqioSi3Gnhzv4OeoIwOJHGa3QF02S4lYhTOQFXVHQN+fpF4mWkp8Haa7Eob584v58w5pSgqbNjTNebtIp1WQ1RRaPVIEbJjgV1tXjyBaMaSd3qthpii8kGdMyP3k+Q+Dc5A1pPHMLqcwTQhxJL+j6mqGlJV9bvEReYQQpwjhMh+v/0o+eIZMzBoBc3uIDtbDm/L52DX6KUsJLmJqqq8uruTAtPY34hHojdUlm/g9X0OAuHYYV9Dcmygqiqt3UHMWdQk6mU0paX1QohThBAPAgeBVsALlACThRBlwIOqqh41Au6FZj2LqvL4sNnH2wecTCvNG5OypNWoY3uzhzPmlKXRSkm6qXcGaO0OUVM0+nJSlz/Cq7s7afeEqCowMqMsjxmleRTl6UedADToNISiCpsaXJw2W76GJjKeYJRITMkJqZJRfQKqqvqOEOJ9YDXxLaIyoB14VVXVnGsOGw1VVj0zy/I40OVnw54uLl1aOeo3s8WopckdwB8em9a9JLd474ATg250VRyqqrKjxcvbtU6iiaig1ROi1RPinQMuCs06FlRZWT61cFT3LrXoeXV3FyumFWOUfQcTFocvnNXpZv0Z9XJYVdUY8XLSf6TPnMwhhODMOaU0u4M0ugLsauth4STrqM7VCAEJSeu5lWNqvJbkCD3BKNuaPVRYRy4n7QlF2bCnq6+keF6lhZUzSmj3hjjY5afO4ac7EOW9gy7CUYXjqkaONEx6LQ5fhC1N3Zw0o+SIfx/J0YmjJ4SSI7UoE3rqhsWo44zZpby8u5O39juZWmIm3zi6P4lWI9jb0SOdwVHKtuZuFFUdscms2R3gXzs6CEUVjDoNZ88tY3ZFXJHUatIxu9yCoqjUdvl4eVcnmxu70Yt8VswcWRG3JE/Py7s6OX5KUdrLWiW5SaMriEGXG6HBhH8Fzq20ML3UTDim8O6B0Vd4FJh07GzxZF1pUDJ2FEXlzf2OEctJPcEI/044gqklZq5eUdPnCPqj0QjmVORz9rz4/v/7DT3s7xi5wMBs0OINRtnRIhVNJypNrkDWNYl6GascxX1CiGMqphVCcMbsUjQC9rT7Rj231qTX4PJHcPnH3qsgyS51Tj9Of3jYfE80pvDizg6CUYVpJWYuWVKJZYSocUGVlZUzigF4aVcHze6RO9WL8nS8vKtD9q1MQGKKSrsnhHkMOaNOb4jOnkhaFqFjjQz+DFwohLhVCHGByHb/9GESi8VQlENlfQVmPUtqCgB498DotGN6f/VGOQrzqOP9gy6MI1RvvLHfQYc3TIFJx/kLRj/QaPnUQuZXmFFUeGF7x4jzL/KNOpy+MLukGu6EozsQGdVWZX/e3tfJ+w0ent/eNu72jDVnsB3YQXwq2YPA00KIRuCPqqoeNa/mLfYP+PUXr2DKnIVMnruYKXMXM3PGInZpBQ2uAI2uAFOKzSNex6TTsLPFw3GTR1dBIsk+3mCUbc3dVA6jTrqzxctHrT1oNYKLF1UMqzKqxGLs3vgG7/17Pe0NBwBQVYVARKHwwlt4Qa/h6hU16IdxPgUmHRt2d7K4uiDr+jSSzOHwhWGE/+7O5np0ej3FFdW0uIPsev05zNVzOGvO+I+IH6szeAUwE48QbKqq9gghzMALwDnjbVy6qN23h2g4yMGdmzm4c3Pf4zXHnQGnXMc7tQauXF494huzwKxjd3sPMWVs3l2SPbY1dYM6tDppuyfEG/vi2kFnzSmlfBin4XF28tDtn8Hd0Zry+claBW8wygd1bgy1r+N1d3HiBZeTX1SadJzVpKOlO0iDM8C0cR5lKMldOr0hhtrtcbQ28uwvf8zujW9w9qduYNV/fY3361wYquZQGGilwDz+8iljdQZe4JxEmSmJCWfvAVvH27B0cuU111FvmYfW3UTj3u007t1B7bYPaN76BroDO9He/AT7O33MqRi+UkivjQuPtXmC1BSNHElIsksskTguGiJxHIzEeHFnBzFVZdEkKwtGKDW2FpdhtljRTtKxctWnWbjyLDQaLYFAALPZTFBfwLM7XXzY4KL7T4/ibm/ipd+v44zLr+Oiz92KRhuPOIQQGLQa3trvkM5gAtHoCmAaUEUWDYd57alfs+GvjxMNhzCaLWh1eppcAZrdQQqmzOPUmba02DMqZyCEmAJ8FogA30msmDXAeaqqnkZc1fSowlxQzMw5M5i/Ij6xs7urnecf/ymmyQtp0+p574CLmWWWkVf8alyaQjqD3KfO4ccdiFBTZEr5/Ad1bryhKBVWA2fMKR30fI/bydPr1nDJjV+npGoyQghu+OFjWIvL0GgOval9Ph8WS7zqaGlNjK1Nbsov+CKT9rzM7o1v8NrfnqBp30dc8637yLMWAVBqMbC9xYPTF6YkzVLaktygyRVMqiTau/kd/v7Ij+hqrgfghHMuYfUNd2AtLuOZD+PR57LJBWmbezCqBLKqqo3Aq4AdeD3xtQG4Ki1WZYHCskqu+dbPuezq6yky6+kORvnXP57lref+OOx5+UbtYesbSTLLewecGIeo6e4ORPpKPM+ZVz5oEdC0/yMeuO1Kdrz9Cs88/MO+xwtLK5IcwUBWzigm36gnXL2M5V/4CTev/S2WwhL2ffgOD9x6Fa0H4yrwGo1AiPiQHcmxTySm0NUTwpTQJNr34Xs8/p0v0NVcT8WUmdy89jd85utrKSgpp9EVoNUTwqjTpDU/OepqIlVVN6qqeo+qqq8nvt4gHikcU2i1Gk6eWUzE3cYbj9/Nc4/+mDf//uSQx+cbddQ7/fil6FhO4w1G2dHqoTgv9ar7vYMuFBXmV+ZTlp98zKYNz/PwHdfg7mhl6vzj+NR//zDlNVJh0Gk4MxFlvHfQReXc47ntwb8yec4inG2N/Pm+/4eSaEEtsxh4q9ZJMCJfS8c6zoQMRe+87dnLTuLUSz/LhZ+7hdsffprZx8UTxKqq8v5BNwAnTC1Ma3PiaFRLX+otIRVCPC2E2JD4eg14M22WZZGZZXlMmTqNkvNuBuAfv7qX15/+XcpjNRqBokK946jR6ZuQ7Gnzog6ROO7whtjX4UMrBCcl+gR6+c9Tv+Ev932TaDjEiRd9ki/d+zsKSyvGdO+Z5RZmlOYRiam8sd9BcUU1X77vSU666Ao+8/V7+yILg05DJKqwvVk2oR3rNLuDKCqEg/HSdCEEn/jStznvM19Epz+0GKl3Bmj3hjDrNSxNlL+ni9HMM7hAPdTh8DBwvqqq56iqejbwmbRalyWEEJw0vRjr8auovPgWAP75xH28tv7XKY83aIXsIs1x3j3owjqEVHVv5/mSGmuSeu3WN17khV//HCEEl33le1xx6xp0hsPbzz9jTil6reBAl58DXT70RhNX3LaGSTPm9h3jcXZSlKfjtb1dxI5AIluS++xt72Hj3x/nods/Q487tfJBPCqIbxsun1o0bHnyeDDWq+uAOUKImUKIXxAfgXlMMrnYRFWBEdPSCznx2m8ihOBfv7mfDX97fNCxRXl6trd4ZBdpjtLVE6LZHcCaooO4wRmI68NoNdimFSU952xrAmDVf93BKR+76oh6AKwmXV938ju1rqQxh6qq8upffsW9N1yMq2EPzp4wB+S8jGMWVVV54sH7ePtvj9LeUEv97i0pjzvY5aezJ94pv7h6dCKaR8JYncGlwH7gb8TLSVeOu0U5ghCCFYkPB8+0M7j8th8ghOCVPz2Kq70l6djeEtNmdzAbpkpGYEeLFwGDPsxVVe2LCpZPKxzUXHb2lTdyy//8mTM/ed242LG4uoACkw53IEJtZ/KHfWfTQcLBAL/9/leIeDp4fW/XuNxTknv86jf/y1vrH0Oj0XL11+9l0crULVpbmrqBeFd7JuYdjPUOu4iXke5VVfW3QPX4m5Q7TC0xU2E1EIgomBaey6fv/Am3/M+fKa4c/GsLBHvaZVVRrqEoKu8fdKYUpdvX4aOzJ4zFoOW4xH6s39uNu/NQE9nU+UvHrStYqxGckJh3sLHe3acvI4TgitvWMGvpCjzOTp5eeys76ttp98jFxbHGrl27+H93fg2Ay7/6PZadtSrlcV09YVq6Q+i1ggVV6Y8KYOzO4CWgDfiCEGI5UD/+JuUO/aODzY3dHHfm6qQ93v4UmXV82NgtVUxzjObuIC5/ZJAyZExReS+xH3vSjOL4bOtImCfv+W8evO0qmmt3pcWeBVVWLAYtTl+Eg/2KDnR6A5/77i8onzyDtrp9vPjgt3h3f2dabJBkB5/Px9VXX00w4GfhGas58aIrhjy2t4hgfmV+xuTNx3qXBmAfsBywMAHmIUwvzaPMYsAfjvFRW0/f45te/QfPP35f38+9KqYjCZNJMsu2pu6UFUR7O3rwBKMU5+mZn5hJ8dIfHqZ26wcA5FmHrueOKSpt3SFa3EFa3EHaPSEcvjA9weiIi4H+0YG9X3QQv2cR/7XmESwFxdRtfYeHfvFzWWZ6DPGnP/2JXbt2UT55Jpd/5XtDRpyhqMLe9vhnzZI0VxD1Z6wf5q8CTUCvVOd04P7xNCjXEEJgm17Eizs72NzgZtEkK57OVtb/4i5i0QhT5i5m2ZkXI4RAALWdvmH1bCSZIxJT2FjvomTAFpGqqmxpjK+8TphSiEYjqN+1lf889RuERsO13/kfiitS74D2hKK4/BFOnlHCjDIL/nCU7kAEbzBGvdNPuyeMdYR31cJJVuz1bjq8YRpcAaaVHJKgKKueyjXf+hmPfesGtrz8NFsOfoOVc4/p3dgJw4033og/FGW7Uk2BdWipm91tXiKKyuQiU0a70cfqDN5TVfXO3h8SMhXHPLPK8ijJ0+P0R9jV5mVxdTWX3vz/+PvDP+KpB+6iZtYCyidPJ9+o5cNGNytnHlMjH45a6hx+ghGFUktyANzkDuLwhcnTa5lbmU8kFOSv938HVVE484rrmbHohEHXUlWVdk8Yk17Df50yjfkp9nF9oShPvF1PY2cAyzCDzvRaDcumFPLuARf2OjdTi81Jq8TZy1Zy5dd+xOQlp/BBU4CT5qhSzfQYQAjBhVdcQ9t7DUP+f6qqyvbmeO4xk1EBjH2byCKEuF8IcZcQ4i7gZ+kwKtcQQvSVHW6q7yamqJy8+iqOO+NiQgE/f/jJHcSiEfJNOhqcAXqC0SxbLAHYVO/GkELHZUtjvEpjSY0VrUlLCRgAACAASURBVEbw7/99gM6mg1ROncWF194y6PhgJEajK8jCaiu3nzs7pSOA+BjV60+eitWoG3G7cEl1AUadhlZPiJbuwYniFedfRlVlBS3dwZTPS44OQqEQ1157Lfv37wfi5aLDfeg2uoK4AxHyjVpmZFi0cKzOYA/wNvHO4y3A5uEPP3aYXWGhyKzHG4pS2+mLV4Dc+n1KqibTcmA3bz77BzRCoEJSYlCSHQLhGNtbPIPCbKcvTL0zgE4jWFxdgKO1kbee+wMajZar7vwJekPyFl80ptDZE+Yq22SuXjGZfNPwwXSBWc/VJ5Rj0seTxENh0Gk4bnJ85Wevd6c8RggBSpRvfm8NGzduHM2vLckxfv7zn7N+/Xquu+46VFVlT3tPUmPjQHoTx4smWdFkWBZ/rM5gMXEJ69cAF9A+/iblJhohWDYl/ubd0hSvGjJZ8rnsK98F4slHV0cLJp2Gbc3d2TRVAuzrSD1nYktTokqjKh+zQUvppCncdM/jXHLTN5g8Z1HSsaqq0tod4qKFFSyfVjTqrZpCk44bT52GRoB7mLGoS2viCpSNriBtQ6z+925Yz/rHf8ENN9xIICCn6h1N7N+/n3vvvReAe+65h0BEod0THHLcqjcYpc7hRyNgUXVmt4hg7M5gB/Bc4vu3ge+Nrzm5zbzKfIw6DR3eMG2e+Kzk+bbTWXr6hZRUVhPo8VBo1rG7zUs4KruRs8n7da5Bbzp/OMaeREVY76oc4nv0p338mkHXaPeEWVxTwFlzy8d8/3KrkRtPm04oqhAaoiLIpNf27Qtvaki9gDjtks9QXD2dvXv3sGbNmjHbIckOqqpy2223EQqF+OxnP8uZZ55Ja3cQIQY3P/ayo8WDCswqtww7nztdjNUZ+IEZQohTgaeAvUdqgM1mm2+z2V5IfH+HzWa71mazffVIr5sO9FpNX1t4774zwCdvuZv/XvcU1TPnx+vVY6qcjZxFXP4wtZ0+Cs3J4fiOFg8xVWV6aR6+5n3UffThsNcoNOv55PHVhx2uTyo08bEllXR4w0OWnB5XU4BGxLcWXf7BeQa90cTlt/0IodHwwAMPsHnzhNmZPar529/+xquvvkpxcTFr164FoMHpZ6jC45ii8lFrPHGcbkG6oRirM3gS0BMfdPMWMHTXxCiw2WxG4ALAYrPZTgNK7Xb774Fim802/kM+x4EliTfvgS4/nkB8CyDPWpikNKgBOeA8i+xs8SA4JA8M8b3/3v3Y46rzefqhNTx8xzVsffP/Bp0fiMQIRVU+t3IKlhR6RmNhxfQSppfm4fSlLiqwGHXMr4qXGX7YmFrscN6SZRx/8dWoqsodd9whGxtzHLfbzTe+8Q0gvj1UXh6PLHe1eVPqYwHs7/ARiCiUWQxUFWSnNH2sr/SFqqquA0jIWp/GkclYXw88AVwOrCIudwHwUeLn94c7uaOj47Bu6nAGCYfC+HxjFwMTwIwSI7WOEPY6BydNPVRZ4ut28dKTD2DML0L5zC0sLxdjno3sdqdOJmabo8UuRVV5eXsrBsDnO7Q9s7czQCCiUJqno+6df9C0byfWkgqmLlye9DqIKSpt3giXLylFE/LQ0TF2NdqBNp051cjj7zvRKXp0Kaqb5pca+Kg1LrO9pNJAnn7wFsHy1Z9j31sv8O677/LrX/+aSy+99IjtyhWONbv+85//4Ha7sdlsrFq1io6ODsIxhf2tLsosOny+wa+BrU3xe80tM+L3D1+AEg6HcTgcmGM9wx43VsY69vJEIUSvxJ4GOI+4QxgzNpvtPOBNu93ut9lsAGXEk9IAQaBqpGtUVIxNV74Xr/BhMHb2jSYcK8un66h1tLCvK8hpcyr62sWdzQfZ9uaLCARzTr+EoG4uM8rGfo/D/b3SzdFgV4PTT5guqvuNtlRVlV0d8ZfW4go9z/zilwB87IY7KC4tS7pWizvI+YsrOGtJ1RHV9ve3qaICLgnreXlXJzUFg4V+LRaYURrkoMPPfmeUk2cO3iaYOdXMyk9/hZd/9SNef/11brzxxiO2K5c4luy68sorWb58OYqiUFUV/xhrcPoxmZxY8wev+rt6wnT0RNBrBYunlIwoP2HoiVBaWkrFOJeejnXs5SbGb+zlF4CHbTbbf4BliZ97fzsr4DiCa6eVCquR6kITkdihfT6AKXMXc8rHPoOixPjP7+5la5OsKso0mxsGy080ugI4/REsBi31r/6JHpeDaQuWcfzZq5OOC4RjmPRazltQMe5NXqfPLqMs30B3IHV1Ua9ExY6W1MUHOq2GBWdcyo8f+R1PPjn05D1JbjBr1izmzJnT93OD0z/k9t7OlszrEKViTGMviTeZhaAvD3Ll4d7Ybrd/2m63n2W3288i3rNwGrA08fRC4MXDvXYm6C0z3drsSdKmv/Bzt2ApKKZ59xb+9syzROSMg4wRisTY1OCmxJIsP7E1UU46w+DhrWfjH6Qfv/mbSR/4qqrS1RPmkqVVaankMOg0XHFCDd5gNOXgmkmF8fkZoagyZL6pJN9IdPJyZMogN3nhhRf4wx/+QCw2uHpsb7sv5esqHFXYndAhWpyFctL+jNUNvUpcwvr6xNcl42WI3W5/GwjabLbrAbfdbn9jvK6dDqaX5lFo0uENRjnQdWiPz2yxct7VXwLgtT88yL42OQEtU+zr9BGJKUkToVz+CPXOAFqNoPGVJ4lFIyw/91KmzFuSdK7bH2FKiTmtA8enl+Zx6qxS2hNlyQM5YUr83lsau1M6DItRh9MXjm8n7d/PQw89lDZbJWMjFApxxx13cOONN/LMM88kPReNKRx0+FI2m+3t6CESU5lUaKQ0P3M6RKnICW2iRHSA3W7/0XhcLxNohGDp5ALe3O9ka2M3s8sP5QZWrvoUbz33exwtdTz46BM8+sM7h7mSZLx494Bz0OprW6KCaF6FhRO/+E3yCwo459M3JR2jKCq+cIzrTpk05oT/WDlvQQWbG90EIjHMAxLFM8ryKDLrcQci7O/0Ma9ysJiZUadhw/Z6vnPVWXg8Hmw2GyeffHJabZaMzOOPP05dXR3z58/nsssuS3qu2R0kmqIBUlVVdiR0iLIdFcAE1SYSIv4B0OoO0pSQIHb5I2Pe0llQZcWgjevL9O8g1ekNXPT52xBCsGP3HilDnAGcvnhvQVG/3oJQVGF3Ystl6eRC8otKuOzL3x000L6jJ4RtWjFTS9KvBZNn0LJ6cRWOFNpFQgiOT2w/DjUboyRPT223yvVfiDu0O++8U5aaZpnu7m5+8pOfAPCjH/0InS55jT1wql0v7Z4QXb64+GH/xWS2GJUzEEL0LmH2AC9wKIm8PU12pZWpxWauW1HJNSuncPmySZwyK14L3tUToXOYBqGBGHSHmtA+bExOFi89/ULu+OWznHb1fw/5YpCMH9ubPYNGW+5q9RKJqZREHRSbUq/4w1EFjRBcsDBz1SzHTymkosDU16fSn3mV+eTptXT1hGlyDZao0GgEQsDpn7yBSZMmsWnTJp599tlMmC0Zgp///Oc4HA5OPfVUVq9ePej5bc3dFKbYItrREl+oLKyypj0iHQ2jjQz+WwgxF/gX0Ai0Jr6OynIZnVZDTaGRhZMKOGlGCasWV/G5lVO5/dxZTCsx0+QK4g+PbjW/dHK8Ca22y59UKaLRaKicNhuzXsPGetcwV5AcKYoSn2Vc1G9ugaKqbGv2oMai1P7+2/zs5ktxtDYOOrezJ8z5CyooNA8ei5kudFoNH19ahSfFMBydVsOSmvgCY2tz6nxTaZ6eTS0h7vx6vLFpzZo1KZOWkvTT3Nzcl7v58Y9/PKgKrTsQoc0TGrR9GYzE2JdYJC7KwLD70TBaZ3AL8Cjw2ICvO9JkV1Yotxr5r1Oncc1JUwhFFVq6gygpEnn9yTfqmJvY293SONg3Fpn1vPTSy9x3/y/SYrMkXjrqDkSS3nB1Dj+eYJTorlfxdDQjNFqKKiYlndcTik86Wzkj8/MnZpVbmF9lpatncHSwuLoArRDUOfwpS1GNei2BSIxTLr6CadOmsXv3bv7yl79kwmzJAB555BECgQCXXXYZJ500WDSh3hkvLhnoJHa1xYUUpxabM7oQGY7ROoPVqqqeo6rq2f2/gI+n07hsIIRg6eRC7jhvNsdPLqSlOzjittHxiSqQXW09BAZEFK72Jp5Zewtr7v4etbW1abN7ImNvcKMfEGZvbfKgRsO43o5/SF5wzVfQapNDdbc/ysWLKrNS2y2EYNXiKsIxZVDlkNmgZU5lfA952xDRgcWg5d16L9/5zncAWLt2LYoiy5gzzd13383999/PD37wg5TP72zxYhrw+lJVlR2J3oLFNbkRFcDom852DvH41vE1J3ewGHVcfnw186qstHmGH1RSajEwrcRMVFHZ3pL85i2dNIWlZ11CNBIZ8gUjOXz84Rib611JvQVdPWGa3UH8W/9FwN3JpJnzWHr6hUnn+RJRwcJJ2XszVhYYOWVWCR3ewaWmvWJlu1pTN6EVmnU0uAKcvfpybr31Vp555hk0muw1LE1UDAYDX/7yl5MazHqJxuI9IwUDBBOb3UG6A1EsBi3TM1C0MFrkq2cYdFoNn15eQ3GefthBJXAoOtje7CE6oCpp9edvQaPV8be//Y09e/akzd6JyO6OADE1/n/Vy7ambpRwAM97TwFw0bW3DPqgdPkjXLCwIum8bHDO3HJ0Gs0gmesKq5FJhUbCMbWvKak/QsR1rzY3efnpT3+a8sNIkj4aGhpwOp3DHtPaHRzU9wKwM6FasDALA2yGQzqDEbAYdVx38lRiqkpPaOhxljVFJsrzDQQiyqA3b0llNYvP/jiqqnLPPfek2+QJg6KovNfgSSon9Ydj7Gn34d30POEeN1PnLWXBSWclnecPxygw6zM+YzYV+SYdFy6sSJk7OK4mvsDY1uRJuVVZZtHzQZ0LX+J1qaoq9fX16TVYAsRLeufNm8dLL7005DG1nT4GalYHwrG+6sJsRqWpkM5gFJRbjXx+5VTc/siQQ2viNeLxN2+qGvFzPv0FtDod69evZ9euXakuIRkjdU4/bn8sSWZ6e3N8ZkFV9WQKy6r6+j364/RHuHBhxaAVW7awTSvCbNAO6keZWZZHvlGLOxChIcV8DJ1WQ0xR+bDRjdPp5Mwzz+S0006jp2d81SwlyWzZsoV//OMfRKNRli5dOvRxTZ5BW0S723tQVJhWYh52/GU2yI13w1HArHILly+rpt0TGrLCaHa5BatRR3cgysGuZBnamslTWHjWJ1BVlfvvvz8TJh/zvFvrxNDv/RSNKX2JuYs+fjnf/M2/mXP8yqRzAuEYVqOOpTXpk50YK0a9lvMWlOMYsBWp0QiWJDpTtzUNUWZq0fPani4s1kKEEHR2drJu3bq02zyR+dGP4kIJN910U58q6UC8wShtniCWfhVuqnpI2DLXogKQzmBMnDijmJUzS/pGXg5Eozk0J3lTw+Do4MwrbuCia7/CT3/607Tbeqzj9kfY0epN2iLa095DIKJQnm+gutCUNHCoF4cvzPkLyrOqDpmK5VOLMOk0g6KDhdXxhqR6ZwBXinnKJr0WXzjG7jYv3//+9wF44IEH8HikJlY6+PDDD/nnP/+J2Wzma1/72pDH1Tv9g5ogW7vjSgd5ei3Tx1l+ejzIrXdEjhMvB6ykKE8/pBTxgiorJr2Gdm+IBmdyaD916hTmrbqeoMacCXOPabY0uhGq2jfNTFVVtjR56H7vKSLv/wm/Z/BgkkAkvqV0/JSiTJs7IqYhogOzXsu8RJnp9iHKTAtMOl7d08WZZ57Fqaeeisvl4tFHH027zROR3qjg5ptvprKycsjjdrZ4Bi04ehPHCybl50TH8UCkMxgjJr2Wz6yYjDcYG1Q1BHGJil71yQ/q3EnRgRACvUbDuwecBINB2tvbM2b3sUQ0pvDmfgfF/cpJ650BHE4Xnvf+xtZ//o6OpoODznP6Ipw/P/eigl6WTy3GqBtcWdS7pbWrzUsoRc4q36il3RMfjvPtb38biEcHMncwvmzatIkXXnhhxKggpqjsavMmzeAORmLsz9HEcS+5+a7IcaaW5HH+gvIht4uW1BRgTkQH9QOig7J8A3//96vMX7Bg2BeUZGj2d/rwheKDaHrZ0tiNZ9PzKCE/s487iRmLTkg6JxiJYdZrOGFq7kUFvZgNWs6ZV07XgOigLN9ATVF8mNKOFNGBEII8g5bX9nRyzjnncNJJJ+FwOPjVr36VKdMnBEIITj75ZL74xS8OOwGttTtIKJJcUrq33UdMUZlcZMqZjuOBSGdwmJw5t4wpxWYcKUoC9VpN31bEB3WupOhAqxEUVU7G4XDy9NNPs2PHjozZfKzw5j4HecZDL91Ob4iGdidee1ywrXeeRH+6esJcsKACY4r5wrnEiumJ6GBABLA84cS2NHWnVNctztOzr9NHmyfEt7/9bYqKijAaszNY/VjlhBNOYMOGDaxZs2bY4/Z3+uLD0hOoqsrO1rgTzxUdolRIZ3CY6LUarrRNJhJTUobuS2qsmPUaOrxh6hwDcgdTajjuvMsBuPfeezNi77FCpzdEbZePon6rqy1NHryb/4kS9DFj8XJmLV2RdI4/HMNq0rN8WnGmzR0zeQYtZ88to2uAxPWUYhMV1ngfS/9Rq73EtyAFb9c6uOCCC9i3bx9f+cpXMmX2hEEIgcEw9BAaVVXZ0uimoF+5c7s3hMMXwaTTMPMwZqJnCukMjoAKq5FLl1bR4Q0NqhzSaw9tSQyMDkx6LcetuhadXs9TTz3F3r17M2r30cy7B51oNYeqNPzhGHuaOvFsHDoqcPoiXLSoImdzBQNZMb0Yg1aT1NMihMCWeD19OMQktLJ8A5vq3XQHolitubsCPdrYvn0711577aii+NbuEG2eEBbjoQj0o4RU9fyq3Ewc93J0vDtymBXTi1lQZaXDO1i/aHG1lTy9ls6eMHWO5L6D6ppqTjg33ncgS01Hh9sf4b0DTsosh1ZmH7UH6NnzLkrAw7QFy5izLLmvoFeZNJ3jLMcbi1HHmXNLB0UHM8ryKMnT0xOKsSeFRIU2MeugVzLd7/ezbt06/vjHP2bE7mOVe++9l/Xr1/O73/1uxGM31rvQaUTfYiUYibGvIyFVnaOJ416kMzhCNBrBJ4+vxqDT9MkC9BKPDlJXFlmNOhZeeA1arZY///nPHDw4uPpFksw7B5yA6NMTCkZi7O4IYFl0Nlfe9RiXfOEbg7qN3f4IqxZX5ky38Wg5cXoJGo1IqlgTQrB8Wjw62NTgRkkhUVFqMfD63i5c/jCvvPIKd955J3fddRehUOpiB8nw7Nmzh6effhqDwcDtt98+7LGBcIyNdS5K+1W5fdTqJZJIHBdbxj7jOBJTCEeVjEyzO7reITlKgVnPlctrcPgjg8L3xdVW8gzx6KB/V7IQgqopU1l+zseYNWuWLDMdAU8gwtu1Dsrzk3MFEUVlakkeK04+jWkLjks+Jxih3GpkUQ7Mlx0rVpOOk2cUD4oO5pRbKDDFu9xTTdDr3Qr79452Vq9ezaJFi2hububJJ5/MiN3HGvfeey+qqvL5z3+empqaYY/9qNVDVFH7FisxRWVronN82ZSRI1NVVQmEY3R4Q7R0B2ntjqubBiIKLd0hWrtD8UqlISRxjpTcEsc4iplfZeX02aW8U+ukpsjU97hOq2H51ELe3O/k3YMuppfm9SkVluTpOeHTt/PV8xcyu/Lo+8DKJO8ecKKoalJUsKWuk3BXCycOkJyA+BvLE4jy+ZXVOb1POxynzCrlrVonsX7D1DUawQlTC/nPXgeb6ruZXW4ZFA2VWw1saermpBklfPOb3+Taa6/lZz/7GatWrcrGr3HUUldXx1//+le0Wi133DH8HC9VVXl9nyNpvOX+Th++cIziPD3TSoZvNI3GFFo9IcryjdimFjOz3MKkQiMleQY0GkE4quAJRugORGlo7aSyYPwrxWRkMI5ctLCSCqsRpy95Nbe4uoACkw6XP5JUCSKEoLiokFf2OORQ82HoCUZ5c7+D8vxDYfbWJg+uLS/T+puv8OHTjww6xxOMUl1kZn5Vbu/TDkeJxcAJU4oGRQcLquLRZpcvPKiPBUAjBIUmHX/f0sIlH/8Ec+fOpb6+nueeey5Tph8T/PKXvyQWi3H11Vczffr0YY9tcgXo8B5KHMeriuKTD5dNLhjksPvjC0Vp9YS4aGEFd5w3m48vm8SSmgLK8o19C0eDTkNZvpFZ5RYWVeUl9diMF9IZjCMGnYbPrJhMMKIkVYJoNYJTZsZHK75f50p6rjhPz8EuH+9s38fXv/51Ojs7M253rvN+XTwq6N33D0VibKl30p2YV1AzZ1HS8b1RwcWLKnNKL/5wOGNOKRFFTcoPaDWHFHLt9e6UC4kCs55Ob4hNDR6+8Y34rOSHH35YzkoeJW63m7///e8IIfr+fsPx/kEXeo2m70O/pTtIZ08Ys17DvMRY3IGoqkpnT5hAROHGU6dxzvyKrEax0hmMM5MKTVyytIr2AeWms8rzqCowEogobGo4pJsjhMBi1PLft9/BQw89xOOPP54Ns3MWXyjK63sdlPZLvm1r9uDavoGYp4OymuksOeW8pHMcvghzKvKZU5G7Nd2jparQxIIq66DhSour4xpYbZ7QoD6WXiqsRl78qJ0LL7mc6dOnU1tbywsvvJAJs496ioqKePnll1m3bt2Ig4N6QlE2N7op7Z/PakyMtawuSDlASVFVWrpDVFiN3HbOLOZWZj+Clc4gDaycUcLMMktSd7IQgtNmxaODLU0evMFDlUdFZj2LLv4cAE8++SQOhyOzBucwG+tchGNKX2I0HFX4sMFF93vrATjtss+j0R4KmWOKSjCisHpJ1bCh+dHE2XPLCEaSK0r0Wg0rEpVFvfmUgRh0GhQFXtnrYO3ataxbt47Vq1dnzO6jnWnTpnHDDTeMeNzOFg+KSt+q3uWPcNDhRyvEkAOU2rpDLJtcyE2nTafkMKqM0oF0BmlAoxFctqyasKImSQdUFZqYU24hpqi8d9DV97gQglkLljDPdhp+v5+HHnooG2bnHIFwjNf2dlGWPyAq2PkmUWczJVWTWXJa8mzjzp4QJ04vorpfEv9oZ2qJmWmlZtyB5NLlxdUFWE06nP4Iu9tSi9JVFBj4sLGbpaedz6WXXopWm9tyHLnA7t27R53DU5R44ri/lPrWpniuYF6lhTzD4L+32x+hKM/Ax4+blFONkLljyTFGZYGRCxeU0z5AzO7kmcVoRFx7v/8g9EKzjqWrrwPgkUcewe0eLME80Xir1kEoqmDsFxVsbnDR/c5fATj7Uzeg1R16E4ajCgLBufOHFhE7GhFCcN78Cnyh5P1+rUawckZcYuODOldKFd3eZPLTm1v6ShLb2tpkwcIQdHZ2csopp7By5Ur8fv+Ix9c7/Th94b5pe4FIrM8xH5einDQSU/CFY1x94mTMKRxFNpHOII2cMquUigIT7n5DSQrMh7ph36519r0phRDMXWpj+qLleDwefvnLX2bF5lyh0xtiw+5OKq2HooItTd0EAwGKpy+kpGoytvM+kXxOT5jz5pdTlJebqpBHwuxyC5UFRjwD5mjMrbBQlm+gJxRj21DzDsx6HL4wb9d5uOuuu5g7dy6vvPJKJsw+6njooYfw+/1MmjSJvLyRB9C8td+BQXtoO3Jni5eoojK1xJyU54J4writO8SqxZVMKc69mSbSGaQRg07Dp06opicUTWpGWz61EKNOQ7M7yIF+jWiFZh2LVsVzBw899BCBQOrE4LGOqqr8Y1sbeq2mL/nWHYiwqb4bjcHM1Xf8kK8/9jy6foJh/nCMfKOOkxNVW8caGo3gY0uq6A5EB83IOGVmPDrY1NA9aFJaL5VWA+/WeVB0ZsLhMGvXrs2I3UcTLperbxH2rW99a8Tjm90BdrR6KU1sY4ajSt8W0bIU8icd3jDzKq2cOqt0HK0eP6QzSDNTS/I4dVZp0paQSa/lpER4/8Y+R1/4LoRgyQkncuLHr+PJvzyF2Zx7q4dMsKvVy552L2WJ6gxVVXljn4OYqjKvMp+aInOSI1BVFUdPmI8tqcx5ieojYU6FhdkV+YOmoU0pNjO5yEQoqrCpoTvluTqtBotBg/m4iykqLubtt9/mjTfeyITZRw0PP/wwXq+3bybESLy6qxOTVtM3bW97s4dARKHSamRKcXLOyheKotNquGJ57jZBSmeQAc5bUEGeQUdPP+2ixdVWKq1GfOEY7x5w9j1u0ms5+7O3sUepSqlbf6wTisR4dmsrxXn6vmqgA11+6hx+XC8+SHH7JhQl+e/SHYhSU2xmSQ4NuU8HQgguXlxJMBpDUQZEB4lKtW0DKtX6U2DS4VX0XHzldQD85Cc/SbvNRwsej4d169YBo4sKmlwBdrZ6KEksWMJRhc2JJrOTZhQnVbIpiorTH+Gq5dU5O9gGpDPICHkGLVccPwmnL9L3JtYIwdnzytAI2NHipcUd7Du+NF9Pg9PPf/Z2Tbgy07dqHXgCUfITCblIYsRlsG4Lnq0v8fwvf0gkdGj7TFFVvKEoly6tytkV13gypdjM8VOK6BzQlVxhNTKnwkJMTa5UG0il1UCh7VLyrVZee+013nnnnXSbfFTw2GOP4Xa7OfXUUzn99NOHPVZVVV7e1YFRdygq2NLUTSiqUF04OCro8IZZPrWI+ZNyW3Ima87AZrNZbTbbepvNdsBmsz2SeOxGm812vc1m+7rNZjumHNW8KisrphfT3m+7qCzf0Dfz4LW9XX15BSEEBSLIlz/3KU5YbpswipNdPSFe3d1FRcGhLaCN9W68wSj+9/4CwJmXX4fRfKiZrMsXZfnUYmbk8NCQ8eb8+RXEVAZVD62ccahSrdGVOt+k02ooKSnixNWfBeCee+5Ju71HA4WFhZSVlY0uKnAH2d3e01fyHIzE+qQnBkYF4aiCKuCChblf4ZbND9yVwHXAYuBcm822AjjDbrf/u97PrAAAIABJREFUFmgHPpVF28YdIQSrF1di1muTtotsUwspMutx+SPY6w+Vk+YXFBLyumhva+XxX/82GyZnFFVV+ee2NrQa+mQnnL4wWxq7CdZvpadhJ3nWQk655Oq+c0KRGFqN4OJFuf9GG09K8w2cPruUjgHRQaFZz4nT47mo1/Z0DbnNWJxnYM45V5JnsWI2mwkGgymPm0jcdNNN7N27l3PPPXfY41RV5eWPOjDpDs0s+LCxm3BMZUqxiZqi5DxfZ0+Y8+aVUZyXG41lw5E11VK73f5y7/c2m20HsArYl3hoJ3Ar8NfhrtHR0XHY989WHf+5M8z8+cNOJhXo+0LMk6dZ+PduN5sa3JQbLFQmjj3j8s/z1P3f4cdr7+Xjl6zO6kzbdP+9trX42HTQQXWBHp8vXjGzYY+bmKIS/iD+Mlj5sc8QU8Hni0s3t3jCnDFZT9DrIjh4EmTWyMRra0GRwoZQCFd3JKlxaV6Jjr1tOpyBKG/t7eDEqYdkDvp/6FcUWbj8nj/x+dPm4PF48HhSl6VmglzqqenpOdS8l8qu5u4QW+s7qbLGX6fBiNInU720ytz32gTwR2JoFZhToBzRZ9VA0vX3yrqEtc1mswINQATofUUGgaqRzq2oOLIV4ZGefziUl6u0BHVsbnAzqTC+tzjLYmFRd4ydrV42toT4VHU5GiFYcc4lvLH+13Q0HuCB/13Pz777tYzb2590/b3auoNsqOtianlBnxrjR61e2rwR1OYduA9ux5xfwNmfvB6TJb4d5PJHmFVl4ZS5+Vn5fxyJTNh06Ql6/rWjnWJL8h71eQv1rN/Uws72AAuqi6gqPPS8xXJoO23uDDMvHfAzb9okKqzZW2hAdt6LAL/97W9xOBx86UtfSvrb9NLfLlVVeX5fPUVWC/mJLaIP9zuIKirTSszMqCxKOtbpCvLZE6cwpXr8CxvS8ffKhX35a4G7gE6gt8vDChyTmVMhBKsWV8Wri/pVfZwyszg+BMcXZWNd3PNrtFouuCY+1Px3jzzA23vbsmJzOglGYvzxgyZMOtHnCFz+CG/si//3a3b/HwCnX/Y5TJa4+mNMUfGFoly2bBK6CZA0HoqVM0ooNOsHVQ9VWI19qqYb9nSlnJcMYDZo0Wk0/PA3z/G9u7+fbnNzDp/Px9133813v/tdXn/99RGPr3cG2NvR0zfJzBeKsj0x37i3VLwXtz/C1BLzkNpEuUhWnYHNZvsE8KzdbvcCLwG9WsQLgRezZliayTNo+bStBlfg0GQ0o17LefPLEcSTpr1T0ZacdgFV0+fS42znhz9/hN1tObQfcoT0Npd19YT6RgLGFJWXPuogqqjMrbBw03fv42Nf+DqnffyavvM6vGFWzihhasnIHaLHMr1Nje5AJKnUFODE6UUUmuO6Rf1zUQMxE+R3d93Mffeu5YMPPki3yTnFo48+SkdHBytWrODiiy8e9thITOHvW1qwGLR9uYIP6uLbmDPL8pIiK0VR8YUVLj1u0lEloZ7NaqIv///2zjtOqvJq/N87fXZnd7b3ZReWKgIClyIgVbBiIWo0QV80JMZXU0zyRvBVkV+iMeXF2GNMLGAhGBTQiKCgIKDIRapUYWEpy/adnS3T7++POzPbZmEpW3m+n89+9s4zz9x77rnlPOU85wBPAytkWd4JTAG2yLL8IyAd6NZZvPuk2Li8VwKnGsQu6pFgZVim1lX9ZG8JlbVedDod02bej8FkxhBws/CrYxwrP3PMlK7A1oJKthypIM1e/yBtzq+gpNpDjMXAhL5JmCxRTJgxC2u0NvZd5/Vj1EtMHdD5hoY6gt4pNkY28VIDzWtocr8kQNNzWZPJ5hDW6BjG3qBNyv/Po/PbVthOhNPpZMGCBQA89thjZ4xwu/G7Mk45XOGJ4KIqN98WOtFJhONDhSh2ehjWw97lGisdOYH8ItA8RdVFxHWXpnGi0kWx0x3O4jUoPYoKt8rh0lpWflvELcMyuHTMFB5+fTUx8UlUubz8c1MB/z2hZ4eP854PhQ4X720rJDXGFJ5IP15RxzfHHEjAyHgXqrsWDPWJQQIBlRKnhztGZGGzdPh0V6fhmoGp7C10UuP2hQOmAWTGWRmYHsO3hU5W7y3h2v6Rx64nzJjFxhVvs3ndGpZ98gU3TT29n3134IUXXqCsrIzRo0dz5ZVXnrZuUZWbVXvrU00GVJXPD5QCMCTL3igEtdsXQOoirqRN6QxzBhctJoOOH47MwqCTwuO+WoTKZOKsRspqvKzdr910MfFaKy/WYkQHvLLhCEVVXdMlsKLWw8KvjmExSuHwEXVeP5/s1bK8De9hZ83fH+cP91zFkb3bw78rrHIzumdCeDxcoBFtNvC9oRmU1Xib5TUYm5eA3WqgrMbD1wWRhxhtcQmMC/YOHnr4kWaRdrsblZWV/PWvfwVg3rx5p+0V+AMqS7edwKyXwi7Pu084Kan2YDPrwzklQBv2LK5yc/2g1C7hStoUYQw6mPgoE3eN7kGVyxdOh2ky6Lj20hSMOomDxTXhaJR+v48tn7zPwQ0r8AcCPP/5YQ6V1Jxu952OEqebl9cfocbjCz8wqqry2f5Sajx+0mLNxBTv4Oie7UhIpOVoWaZKnR4yglnkukvSmgvJgPQYhvawU9zkRW4y6Lj6khRtMVqJi4PFkfMeTLz1R1htsRzdtZlHX/pXiyEtugOLFy+msrKS8ePHM2nSpNPW3XKknCOltSQ0mDT+Kl8LH3NF78RGbr3lNV5yk6IYkds1gyUKY9AJyE2M4uYh6ZyqcodbdgnRJqb0TwZgw3flHCiqJn+XwpIFj/Cff/4fVtVDlEnPKxuO8E1B5/HTPh2FDhd/W5+P2xcID4upqsoX35VzuLQWo17iyr4JfPzGMwBMueNeLFHR1Lh9IMEPR2Z3qmQgnQlJkpg+KB2DXkedp3Hk0uQYM1f01iJlrt1f2iikeoioGDsTb70HgM+W/J1Fmwsa5eruTtx7773861//4sknnzxtvco6Hx/uKiIlxhRugGw8VI7Hr7mS9kqqnxPw+gO4fQG+NzSzy4ZFEU9WJ2Fkz3jG9EqgyOkNhyjunRLNqNw4VLQJ5UD6QHIvGUqt08GGFW9iMxtIijayWDnOmn3FzTxKOhMF5bX8bX0+QLiVpaoqmw5XsPNEFToJrh6YwqGvP+HUkYPEpaRz+bXfx+cPUF7r5QcjssKhggWRibEYmHFZBqXVnmbDRZdmxJAbb8brV1m1pzhiIpxxN8xk8vd/wt2PLKCgvI5/f3OiRbfUrowkSdx4443IstxiHVVV+Xh/BUiEhzKPV9RxoLgGvU5iQp/ERj3Uoio3Vw1MDc8rdEWEMegkSJLEdYPSyI4zU+ys9/wYkRvPiJyQQSjl0htmA7Du369R46jAbNSTHmvm42+LefebE1orupNxsLial784gtmga5R4ZvORSrYdc6CTtEnQzBgDqxZqKT+vmvkAeqORQoebqy5JoV9axycM7woMzopleE5cs3F/SZIYmxtDrMVASbWHjYfKm/3WZLFyzaxfEBOXSIbdzLZjDj7ZU9xtsqKdOnWKgwcPnrkisOlQOQdL6kix1bs8f35AW/syIieO2AbRRytqPaTZLYzN65rDQyGEMehEmAw6bhmcRKLNREkDV8GRuXHIQYOwW80ha9BoXLXVfPKW5oxl0OvIirew47iDBWu+Y/8pZ6d4gF1ePx/uLOSVDUeIMeuJaeABtOVIBcrRSiQ0z4ueSVFsXvkuFUUnSO2Rx9BJ13PS4WZgRgwT+yZ33El0MSRJ4sYh6SRGmymvaexOajLouCo4f7DrpPO0a1b8Pi/GqhOs2V/Cp/tKOsX9dL7MmzePoUOH8sYbb5y23u4TVSzfWUiKrT6M+lf5FVTWeYmzGhs5MPiCaSxvGZYZnmDuqnRt6bshUSY994zJIdpspKxaG9uVJIlRuXHIPeyoQGDETCRJx5cfLaH4uDb0opMk0u0WDDqJf246yvvbC6n1RM561R7kl9bwzNpDbDhUTobdEnZ5VFUV5Wglm4OrrKcOSKZ3sra2IiU7j4xe/bnqrp9T6PTRJ8XG7XJWlx2D7SgsRj0zR2Xh8avNMp+lxpoZF5w/WLOvlCNlzdesVBSf5M8/uYF/PnovyVZYvaeYFTtPdekho507d7Jw4UIAxowZ02K9o2W1vL3lOEnRJozBdJYHiqrDPdjJ/ZIa3Y+nqtxM6pPcKdNYni3CGHRC7FYjs8fmYNBLVNQ2MAg94xnew44xOZfowVMJ+H18sezNRr+1mQ1k2i0oRyp4es13fHuyql2T5Li8fj7YWcjf1ufj8wfIjLOEHx6fP8Cn+0rD8fan9Euib2r9OoI+Q0fz82eXEH/JWPqn2bhzVHa3zlzWlqTGWrh1WAbFTk+zuaTBmbEMy9YaFiu/LW6USwPAnpSGOSoKR+kpvl75LzLjLGz8rox/f3OiSyZcUlWVOXPmoKoq9957L3369IlYr8Tp5rUvC7CZ9eFk9SVOd9i9e1xeIhlx9XGeSqs9ZMZZmdy/e/RchTHopCTaTPx4XA4BVQ0nQZckidHBiea4cTOJv/Je/KP+q5kboE4nkR5nQQcs+uoYf1x1kA3flTYKnX2hKapy89HuIp5ceYBNh8pJt1sajatW1XlZuq2Q/UXVGHQSV12SzIB0bR4g4Ndar6qqcqrKy4D0GOE5dAEYkmVnTK8ECiOsG7i8VzyXpMfgD6h8uKuI0ur6OjqdjmtnPQjA2sV/x11TRWa8hW8KKnlny/Eu52W0atUq1q5dS1xcHA8//HDEOk6Xj9c2FSBBeDjT5Qvw0bdaaJT+aTYGZdbPW9V6/ARUlR+MzOo292n3OItuSmqshR+NzcXlU3E0MAjDeti57YpLyBhzE8U1fhYrJyJ2920WA5nxFkx6iQ93FfGHjw+wbHshBeW1F6TL7/EF2FNYxcvr83l6zXds+K4Uu9VARoPeAMCx8jqWbD1JSbWHWIuBW4Zl0Celvkew9Ln5vPPnORw4eoIB6TH8YIQwBBeCkFNCut1CaZNwFJIkMbFvInnJUXj8AVbsKGrkctpPHkfe4BHUOh2sWvQcOkkiM87CnkInr315NNxA6ez4fD7mzJkDwJw5c0hMbJ6M3uX1s2hzAVVub9jTLRBQ+fyQA6fLR0qMiYkNvIf8AZWyGg+3y1kk2bqu91BTxBPXycmOt/LT8bmA1MjLKN1u4fvDM8lJsFLrKOPfKz/XFm5FaP1bTXoy4ywkRRtRjlbw0rp8nly5nw92FpJfWhPRzbAlKmo9bCuo5I0vjzL/P/tY+OUxipxuMuxm0u2WRi9xrz/A5vwKVuw8hcsXICfBym3DM8IZogCO7t3B16uWsn3dSnrFStwxovu0tDoDJoOOmSOzMegkqpr2ICWJaQNSyI63UOv1s3znqUbDkjf+9GF0Oj1f/udfnDi0F0mSyLCbOVZex7OfHaagC8TIeu2119i3bx89e/bkvvvua/Z9rcfP618WcLy83nMIYNPhcgqrvFiNOq4dmIohODmsqiqFDjcT+iQxMKPrRCRtDSLASxcgM87KAxN7sXBzAYUOF2mxZiRJwmrSI8c62fSPe1ENVnalvsyBomouy7YzNNve7KVq0OvCftBuX4DNhyvYeKgck0FHQpSRWIuRGIuBuCgjNrMelzdArcdPrcdPUXklbhzh3LtRRh3JNlPEyV1VVfmupIZNh8pxurUhIDknjpG5ceE4RKANDy194fcA3Drrp/zy5iu6VJTHrkKizcSPxuby9Ko9VLt8jeI6aZniUlm+4xRFTjfvbj3JVZckk5MYRXrPvoyZfgcblr/J+y/+nvv/8iaSJJEaa6aqzsuL6/K5cUgao3smdNpV4ZMnT+amm27i1ltvbZYcqsbt442vCjheUUea3Rw+h28KKtl+vAop6PLcUF+lNV5yEq1dMvbQmRDGoIsQF2XkJ+NyefebE+w6XhUeikntkUdqVi4nvtuDfteHeIfdwpajlXx70smI3DgGpNnCrZqGmA06UoPRQn3B1ZOFjjoKylW8gQA+v4peJ2l/koTb7SYuxkhGg4cmEiVON198V85JhzYpmRRt4oo+Cc3SAQKsXbGEwkN7SE3P4MU/zReGoA3JiLNw+9Bklu5xotNJRJnqJ+ZNBh03Dknj030lHC6t5YNdRVzeM55hPexMu/MB8ndvZez1dzTaX6zViMWo571thRSU13HjkPRwPorORF5eHosXL25WXu3W5ggaNq5UVeWr/Aq2Fmj5jMfkxDSaMK52+5CA2+WsLu9GGglhDLoQFqOeH4zIZlV0MZ8dKCEp2kSUSc/0H/+Wvz00i+Ofv80Pr57O3rpYipxu1h0s46v8CvqkRNM/1UZqbOQXuUGvw6AHaPlhrpG8jV4gDfH6A+SX1rK/qJqC8jpUwGLQMTo4Salrckx/QOXoySLWv6MtMFvwlz9HzDIluLBkx5m5a5Sd174sQC+ZGnlqmQw6rhmYwpajlXx9pJIv8ysorfYwuX8Sv3ju3Yj3jcmgrW/ZfsxBQXkdtw3PJCexc4RtLi4uJjExEb2++T1b7fLx6qajFDvdpAcbRKqqsv5gGbtOOpGAKwckk2WrP+cat48ql4/ZY3IaRSntTnQ/89bNCSWAnzW6B3VeP6ccLnoOkhkx9WZ8Xg/rX32CGZelcvXAFJJtJty+ALtPOvn3tkLe+voEW45UcKKy7rxdBL3+AAXldXy6r4RXNxWwem8JR8vrkCQYkhnLzFFZXJoR28wQVNZ6OelwceA//6DO6WDixInMmDHjvGQRtJ7+6bF8X86iuNrTzCtIkiRG5sZrQRL1EgdLanj3m5ON5qpqnY3jYOkkiYw4C26fnxfXHWbFzsJmsZHaG7/fz4wZM5g0aRJHjx5t9F1RlZuXNxyhuNpdH5I6oPLpvlJ2ndTyE1xzaQr9Grg8V7t9OFx+fjQmh7wGjg/dDdEz6IJIksTAjFh6JESxYmchO445mDLr1+xV1pP/7TdsXvkuY66/nd7J0ZRVe9hXVM3+U9VU1nnDi70ktGB4qbFmkmwmok16okx6ok2aj7VekvD4A3j9Kl5/AEe1lxpHFcVON0VVHsprPDT0R0qNMdMv1UaflOiwj3ZDPL4AxU4PaXYzd47K5v8+s/GFycTTTz/dacebuyvDesTh8vpZvqOQ+ChjoxwIAL2Sorl1mJH/7C6ivMbLu9+c5NJ0G9Wb3mLT8jf5+V/fIS23sa++3WrEZjbw5aFydp+o4tbhmfROju6Qa/v888+jKAqZmZnEx2uJZ0KLHZdtL8RkkEgN5gJxef2s2VdKflktBp3EdZemkp1QP6RZ7fbhdPmZPTaHvOTu3XsVxqALE2Mx8IMRWQzOiGXp9pNM+K/f8p+/PsTHbzzDsMnTsURFk2gzMdaWwOU94ymoqCO/tJZip5vSag9lNdrfuSABSTYTPROj6JdqaxRzqCH+YEIaSYIbBqcxsmc8Rr2O5557jrlz55KRkXEeGhCcK2PyEkmymXjr6+PUeT2NPLxAayjcLmey5Ugl24872F1YTdXhQrzuOt57/nf89I+voWsyBKPXab0Ep8vHKxuOMDjDzsR+SWTGWdrNKBw+fJjHH38c0IxCbGwsNW4fy3ZojaaUmPrhsaPltazdp4VONxt0XD8olXR7/RxBjcdPwO9n9rgceiV1b0MAwhh0eSRJYlCWnZzEKL7oEUfx4T30HnklqsHSqJ5OJ5GbGEVucEzX6w9Q4vRQ5HRTUeuh1hOg1uMLew+pKhiDCT1MBh16SSUh2kxKjDncmzjdJJrbF6Cs2oMkSYzIjWNi3yTio0x4vV4I/k4Ygo6lb2oMP5+Ux6LNxzhZ6SLNbm40rGfU6xiTl0DfVBufHyjFP+aHOPZsIP/brXz49j+44c57I+43xmIg2qznQHE1u05W0S81msn9U8hJsLapUVBVlQceeIC6ujpuu+02rr76ag6V1LBk63GqXX6y4jWj5PEF2HS4nN3BZPapsWau7J9MfFTD4HNear0BfjYhN/zMdHeEMegmxFqNXDc4nYmLXmD7MQefHSilrMJFtFlHrMXYzAXUqNeREWdp5C0RIhSUrOGDW1NTc8ZJXlVVcbp9OOt8WIx6pl2SgpwTH17R+eWXX3LPPffwyiuvMG7cuPM9ZcEFINFm4qfjc1m2vZCtBZWkxZqbuSQn2Ux8b2g6ewpjWD39QU7+ax5fLH6R6sQBXHnFqIjpV3WSRHKMCVVVKSiv46V1h8lJjGJCnyT6pNjaZC3JwoULWbt2LYmJifz0ofm8uC6fgvJaYi0G0uxmVFXlWHkdnx0opcrlQyfByNx4hmXbw55s/oDKqSo3yTYzN8upF40hAGEMuh3RZgNjeycysmc8/3jnPbzWPIqqVQIBzVU0LsqI+QwP4tm03vwBbXV0nTeABGTGW7n6klQuzYht5K1SVVXF3XffzZEjR1i5cqUwBp0Ii1HPbcMz6ZFg5cPdRRgkiaQGETshNE8VQ94Pb+b1Uzs4tO49dr31BKWGv9I7Q0tFmhbBW02SJBJtmlEorfawaPMxjHodw7PtDO0Rd8ECvBUVFfHbhx4CYOJdv2bF/lqizdpiy4BaH2wutE4mMdrE1AHJjYbHqt0+Kmq9jO+dyNQBKTgqyi6IbF0FYQy6Kc88vYBHHnmEa665hrcXL+FYpZu9p5zsOO7QJn9VCRUVq1GH1ajHqNfWFEQyBKqqElBVPL4AdV4/dR4/fhV0kvaw902JZkiWnV5J0Y3iETXkwQcf5MiRIwwdOpR58+a19ekLzhKdTmJMXiJ9U20s367FkEqKMWFtsnbAYtTzowf/l6cP76Dk2CEcGxZxePKPOVxai91qoF+qjX6pNuxN7gNJkoiPMhIfZcTrD6AUVPLVkQpiLAYyolQG1xlJjbWQEtO8Z9ISHl+AQoeLI2W1fHvCyWXT7+HEgR0MHH8tMRYjLq+f7cer2HncEV78aDHquCxLW5Spb9AbKHZ6iDLp+Mm4XHp3Y4+h0yGMQTfllltuYcGCBaxcuZI//uEJHn/8cfqm2pg+KA2Hy0tZtTZfUFBey/EKF5V1PtxBV0NJAgltEU7IY8jj9pFoCJBht5KdYCU11kxClIlkm+mMkUWXLFnCW2+9hdVq5fXXX8dk6p5+2t2BJJuZu8fksP24g+U7CnHUeUm2mRsNMxrNFmbO/TMfvbqAa++9n3yXlf1F1TjqfHx9RFunkBZrJjcxiqw47QXfcEGhscFKeJfXz96iag6Un0SSQEVrtdvMBqxGvfZn0hor1W4fNW4/Trcv7PcfishqM+u59vZ7cHn95JfVcqikjOOVdYRCcMVZjVyWHUv/1PpFmD5/gNJqL35VRc6J4+qBqdjMF+8r8eI9825Obm4uixYtYvr06Tz11FNcdtll3HTTTeh0EvFRJuKjTPROsTE2rz5wl6qqeP0qHr+2AtmgkzDqJQx6HWWlJaSknP0S/IKCAn72s58B8Kc//Yl+/fpdsHMUtA06ncSwHnH0To5m1Z5ivimoRJJo5DSQ0bMfs3/3sraNFgX1eEUd+4tqOFxSw6kqN6eC0VKNeokMu4XMOAsJ0SbirEZiLQZ0OgmLUU9ClJHoaG3uSlVVXN4A5TUe/AEVv6riD6ioquatZNBJGPTa/6RoI9s3fYYUn4U3OoVCh4uTla5wA0YCsuMtDM60k5tYP3nt9QcodXpAgtE9ExiblyhSqiKMQbdmypQpPPHEE8ydO5fZs2fTr18/BgwY0GJ9SZIwGaQLNrmnqiqzZs3C4XBw/fXXM3v27AuyX0H7EGs1cuvwTKb0T+ar/Ao2HirDH1BJjDY2Cj0R8Pv5YtkiLr/u+/QYkIynTyIF5XUcr6wL9jq9HC2v42h5Xfg3OgliLUbsVgMGSSXG6sJi1GM16sIt94YDlu7gEGXI263G46fw8F6Ov/E/SHoD6fc8hyFWy+KWHW+ld3I0PROjwmte/AGVyloPLo8fg0HHhL5JjO6V0Gw462JGGINuzi9/+Uu2bdvGkiVLuOWWW9i4cSNxcXHtcmxJkpg7dy5z5szhpZdeEovLuigJ0SauvTSV8X0S2Xq0knUHSymv9YKqxcxa/c+n2Ljibb7bsZlZjz2LyWCkd0o0vVM077Nql4/jlXWcqnJTWeulss5LtdtPZZ22reFqWYAI+GsdFL77e1Sfm5TLJjFyYG9SYixkxVvChsrr13oYdR4/Op1E/7QY5Jw48pKiRdKkCAhj0M2RJIm//e1v7Nu3j5qaGsrKytrNGABMnTqVyZMnR4wRI+ha2MwGJvRNYlzvRC2kSJGTbwoc5I2fwda1/2HflvW8/eeH+cFvn2p0vW0WA/3TYuifVp8cxusP4Kjz4XT5cNTUEpC0Cd86rx9fw1wbqjaPYDLoiDLqsJr0WPQqq/5vPv6qYrL7DeK+R59EZzDh8vqpcvkor/GiAlajnl5JNgZnxtAnNabF2FoCDWEMLgKioqJYunQpNpstvDy/rVBVlV/96ldcd911XHnllQDCEHQz9DqJ7Hgr2fFWJvdLpmxMD8bnLuaBu25h5/qPkCw2Jt/9W1Qk9DowG/RYDMHFi8GJZKNeR5LNRJLNRI1VPeMalkBAm8uqra3l3b88xLFvtxBlT+TqX/yZkjoVneQhJcbMJemx9EyKJt1uJjHaJHqjZ4EwBhcJ2dnZ4W1VVXnhhReYOXPmBe0lqKrKww8/zEsvvcQ777zD/v37sdvtF2z/gs6HJEkk2czMvH4SWe8t5cYbb2TH6iWMGdCDmff/hlKnm2Knm2Knh5JqbVJYIuixJoGqgsvlxeJzUT9LUD8FHPJo0+skYsx63n3yZ3y3cwsx9jie/scixl4+gvgoU3hCWnDuCGNwEfLMM88wZ84cXn/9dZYvX05mZuZ577Ouri5sCAwGAwsXLhSG4CJj4sSJLFq0iNtvv52XnvkLk8aO5IYbbgh/H/ZW8wXw+AN4fAHNs6esjKTERJAkJLTVyw00hvk/AAANVUlEQVRDoRiD3kOSJJFbeT+PPfYYK1asoH///h13st0QEcL6IuTmm2+mX79+7N69mwkTJrBv377z2t+uXbsYO3Zs2BC8+uqrTJs27QJJK+hK3HDDDbz88ssMGTKk2T2geavpsFkMJESbSLNbyE6IItNuJjshiux4K1nxVjLiLCTHmImLMhJl0hPwecPDPbfffjvbt28XhqANEMbgIiQnJ4e1a9cyevRojh8/zhVXXMETTzyBw+E4630tXryYsWPHsmfPHvr06cO6deu47bbb2kBqQVfhzjvv5IsvvsBi0dYOOBwOXnvttXDMq9ZSUVHB/Pnz6devH1u3bg2XR0VdPPGC2pOLdpho9YEKnnt1b9i1TYJG8fmjTHpMeglHnY90u4VfTclj+pD0Rvv4YEchC9YcotDhwmLU4fIFtMUxEozMjeNouYuTDhd6CfzB8ob/M1rYb3uQmJjIRx99xOzZs3nvvff43e9+x/PPP8/GjRvp1avXaX8bCATQ6bR2xODBg9HpdMyePZsJMx/kN+sKKVzxaat0FqlOa79vqNdIeoy0H6DZb+OsBjx+ldpgQpY4q4FHrukX3lfD/ditRmrcXrzBnDCSBLcPz+Tx6S2v3Wgob2iMPITVqMNs0FFZ52t0LhP6JLLuYNkZz/FC0JK+H/9gL0u2ngjfr7e1cJ4t/b7hKvNf/OIXLF68mPfff58HHniACRMmNMtH3BCHw8Hzzz/Ps88+G26gvP322wwfPvyCnvuZWH2ggn++eSB8bqHr0tK9GaIlnUS6d0OEdDy8R1yr7u+2QDpba93WyLL8a6AYsCuK8nxL9ebMmaPOnz//nI7xwY5C5i7bgzfQ+nO3GHX8fvqARi+JRz7Yi8t7fhnDmu63uLj4nFb6ng/r169n/vz5eDwe1q9fjyRJOBwOXn/9ddLT00lPT0en03H48GE+/PBDdu7cye7du8NeQkePHmVnpamZPlqjs4Z1zuX71tYz6rXJSF8rLpdBB0/dNBCgVdf4DrnxizJ0DS/UPRKiqT7Plqb3Vkv6HpoVy5f5lc1+3/Q8z3S9Qrz33nvcf//9VFRUABATE8O0adO47rrruP7663G5XKSkpLBu3TqWLVvG4sWLw3UnTZrEo48+ypgxY87pnM+VD3YU8r8r9uD2tfyOiHSuLenk5iHpvL+j8Iz3QtNGaaRjne87Yt68eTz11FPNZts71TCRLMvjgERFURYB8bIsj2qL4yxYc+isDAGAyxtgwZpDjfZxIR7ypvvtCMaPH8+nn37KihUrwmOz+fn5PPTQQ9x1111MnTqVKVOm8OMf/5jly5eTn5/Prl27wr/PycmJqI/W6KxhnXP5vrX1vP7WGQIAX0DbR2uv8ZKtJyKWX6h7JMSFvlda0nckQwDNz7M11xxgxowZbN26lblz5zJ48GCcTidLly7lnnvuYfXq1eF627Zt46WXXqKiooKxY8eyevVqVq5c2e6GALRzO50hgMjn2pJOlmw90ap7oaUjtsd7orMNE10L7A1u7wl+3txS5eLi4nM6SKHj7FY7Nvxd6Jjnuo8z7beyMvKD2F6E5PB6vdx9990UFRVRVFTEqVOn6NWrF9OmTWPq1KlkZGQ00n9L+miNzkJ1zvX7s63XGs5mH3618b0YuoYX8h5pKNe53vdN762zla/pebbmmocwGAzcd9993HfffRw7doxPPvmE1atXU1hYGJZr0KBBzJ07l6FDhzJq1CgkSTrncz1fWqubpufa0u/8F2AAJnSstnpHdDZjkARUBLddQNrpKp9rVyndbuHkOTyo6XZL+Jjnuo8z7RfO/bwuJCkpKYwYMSL8+Uxd05b00Rqdheqc6/dnW681hNIftmY/eqn5NbtQckSS63zuj4a/PVv5mp5na655SzIMHz6cOXPmAPX31pQpU5gyZUqr5WlLWqubpufa0u+azhGcq0yhY7XFO6JTDRMBJUDIVSAGaJPsEr+akofxLBeoWIy68CRkaB8W4/mrr+l+uyqR9NEanTWscy7ft7ae5qveunMx6LR9tPYa3zY88jqNC3WPhLjQ90pL+r68Z+SFiE3PszXXvKvyqyl5mA2nf0dEOteWdHLb8MxW3QstHbE99NrZegYfAdcAS4BLgI/b4iDTh6TjcFbx3IZT5+xNFNruqt5EF5qm+miNzprWOZvvT+dt0dJ+Iv32TN5EDfdztt5ETeXtbN5Ep9N3a7yJWnPNuyqhd8Q/vy45K2+i0+kkkqdQCOFNFAFZlh8BTgDxiqIsaKne+XgTQcd47bQGIdfZ0Rnl6owygZDrbOmucrXkTdTZegYoivL7jpZBIBAILjY625yBQCAQCDoAYQwEAoFAIIyBQCAQCIQxEAgEAgHCGAgEAoEAYQwEAoFAgDAGAoFAIKATrjM4G+bNm9fRIggEAkG3oNOtQBYIBAJB+yOGiQQCgUAgjIFAIBAIhDEQCAQCAcIYCAQCgQBhDAQCgUCAMAYCgUAgoIuvMzgdsiz/GigG7IqiPN+gvC/wfaAW+EBRlAORyjpArjuAXwKxwJ2KoijB8ueAW4HtiqJc3d5yBb9bBoxG082PO4m+VgH90RLUqYqi9IwkaxvKNR6YpyjKlCblY4CxaA2t1xRFKY5U1gFy/RK4C01ftyiKkh8s72h96YCvgB7Ai4qi/L+O1pcsyxKwC7AFiwoURRkfSdY2kikGeBUYDnysKMp/N/guBXgAOIX2TtgUqexcjtstewayLI8DEhVFWQTEy7I8qsHXzwBPA88DT52mrN3kCt58tYqijAL+AswPlmcC3yiKktbGhqBFfcmyPAJ4KShD6GXR0fqKAX6lKEoOmkF45zSytgmKoqwHrBG++gPaNXyH4HVsoazd5JJlOQ7tJTEMLbXsb4LlnUFfM4C7gjKEXq4dqi8gG5imKEouMA5YdhpZ24LRwCzgUmBK8DqFeBJ4U1GUF4G5wXdHpLKzprv2DK4F9ga39wQ/b5Zl2QrkKYpSDSDLcs/gi6VpmUFRFF97yaUoigosD5ZvAS4Pbk8GHpVl+XvALEVRSttAphblCn6eBPxMluW1wH1oLcuO1pcT+DZYPg1YHUlWRVFq20Cmhngafgj2mHzB61kgy/IVkcraWKZmcimKUgl8Hvy4BRgc3O5QfQUZCzwny/JbwG+B3nS8vgoafLwZeD+SrIqiBNpCGEVRPglty7K8G63FH2Ia0NBw57ZQln+2x+2WPQMgCagIbruAtOB2PFDVoJ4PbVimaVlyO8vVkCuBBQDBFnEesCZU1t5yKYryJ6AnUArMIbIOO1JfVwBftCBre9NQXtB0FamsI7kc+Dt0Cn2hKMqDaPd4FlpruLPpq2doSC2CrG1KsKFaoCjKsQbFxqChhPpnIlLZWdNdjUEJEBXcjgHKgttlgKVBvSigOkJZZTvLBYAsy72Bo4qi7AmVKYqiKoryNGBqI5nOKFew1f8Q2osjkg47Sl8GwK8oir8FWdubhvICuFso6xCCww2rG47Bd7C+QjLUAj8HLqNz6SsZKGpY1kTWtuZO4LEmZdUNtkPPRKSys6a7GoOPqO8KXwKskmXZriiKGzgqy3KULMsW4JiiKI4IZXXtKReALMupwBBFUZbKsmyTZTk6NPYny7IJrXvfVpxOrtD4YwywoQUdtru+gkwCPgt9aCprG8nUDFmW9bIsxyiKcpCgoZRluRfweaSy9pYruN0HiFMU5TNZllNkWZY6Wl/B7ZAMicCazqKvIDdSP3zbTNY2luUmYJmiKE5ZllNlWQ71kD4PNhoBzEHnjUhlZ023DVQny/IjwAm0buYaYI6iKHfIsnwpcAtai2O5oih7IpW1p1xongBr0IZcACRABpYANWjj9wsVRalpT7mC+toIbAv+va4oir+j9aUoyh3B755E8wTxBj83k7UN5RqEZqyuQWtVj1MU5SFZliejXTsL8IqiKIWRytpTLrRhoY8AJ9q9VagoyvUdrS+0yeGvgI+BbYqihBwBOlRfiqI8FPzuL4qihCbboyLJ2kYy/TfwP2gtfBPwJnB90KMpA/gZ2jzCNkVR1kcqO5fjdltjIBAIBILW012HiQQCgUBwFghjIBAIBAJhDAQCgUAgjIFAIBAIEMZAIBAIBAhjIBCcN5IkfSJJkniWBF2a7hqbSCC44EiS9CxaqIR+QArauodJwFWqqrZJnBqBoL0Q6wwEglYiSdIAVVX3SpI0C+ivquqcUFlHyyYQnC+iZyAQtJIWXvpGSZI2ABOBP6JFwLQCA4EX0PI+fKuq6u8lSRqJFpVzKrBcVdVlEfYnEHQIYpxTIDgPVFXdGfzvQwu3Xaiq6i/RgvcdVlX1DuCqYPXfAOXAWjRjIRB0GkTPQCA4f3wN/ofCe9c02A7F/LlEVdWPAcSEs6CzIW5IgaD9UCVJmhHcvrZDJREImiB6BgLBWSBJUhwwBsiTJCkLLapkniRJecBQtBd+FtAHGC1JkhXoJUlSL7TotAslSXoQ+K+OOQOBIDLCm0ggEAgEYphIIBAIBMIYCAQCgQBhDAQCgUCAMAYCgUAgQBgDgUAgECCMgUAgEAgQxkAgEAgEwP8HmrlMTk57pGoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Get the average predicted intensity function, and the intensity confidence region\n",
    "model.eval()\n",
    "with torch.no_grad():\n",
    "    function_dist = model(quadrature_times)\n",
    "    intensity_samples = function_dist(torch.Size([1000])).exp() * model.mean_intensity\n",
    "    lower, mean, upper = percentiles_from_samples(intensity_samples)\n",
    "\n",
    "# Plot the predicted intensity function\n",
    "fig, ax = plt.subplots(1, 1)\n",
    "line, = ax.plot(quadrature_times, mean, label=r\"True $\\lambda$\")\n",
    "ax.fill_between(quadrature_times, lower, upper, color=line.get_color(), alpha=0.5)\n",
    "ax.plot(quadrature_times, true_intensity_function(quadrature_times), \"--\", color=\"k\", label=r\"Pred. $\\lambda$\")\n",
    "ax.legend(loc=\"best\")\n",
    "ax.set_xlabel(\"Time\")\n",
    "ax.set_ylabel(\"Intensity ($\\lambda$)\")\n",
    "ax.scatter(arrival_times, torch.zeros_like(arrival_times), label=r\"Observed Arrivals\")\n",
    "None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
